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1.0  INTRODUCTION 


Effective  transonic  operation  is  required  of  modern-day  transport  and 
fighter  aircraft.  This  consideration  has  stimulated  development  of  numerical 
methods  to  compute  transonic  flowfields  about  increasingly  complex 
geometries.  In  order  to  be  useful  for  design  and  performance  analysis,  such 
methods  must  contain  accurate,  reliable  schemes  for  solving  the  nonlinear 
equations  governing  the  flowfields  as  well  as  the  capability  to  handle 
realistic  configurations. 

Various  finite-difference  methods  have  been  developed  to  calculate  three- 
dimensional,  transonic  potential  flowfields  about  isolated  wings  and  wing/body 
configurations.  To  date,  these  methods  have  used  formulations  that  require 
alignment  of  computing  coordinates  with  appropriate  geometric  surfaces  in 
order  to  apply  surface  boundary  conditions.  The  codes  of  Ballhaus  and 
Bailey*-^  and  Boppe,^-^  which  solve  various  forms  of  the  transonic  small- 
disturbance  potential  equation,  use  mean-plane  and  nearest-point 
approximations  to  represent  wing  and  fuselage  geometries,  respectively;  the 
approximating  surface  is  chosen  to  coincide  with  a  cartesian  grid  and  is  the 
location  where  a  linearized  f low-tangency  boundary  conditions  is  enforced. 

The  codes  of  Jameson  and  Caughey7-^  are  based  on  the  full  potential  equation 
and  require  application  of  an  exact  surface  condition  using  analytic  and 
numerical  mappings  which  generate  computational  grids  that  conform  to  the 
entire  surface  geometry.  Although  the  small-disturbance  approach  has  so  far 
produced  the  most  extensive  geometry-handling  capability , comparative 
calculations  indicate  that  methods  based  on  the  full  potential  equation 
provide  greater  accuracy  in  matching  experimental  data.*® 

This  report  presents  an  alternative  approach  to  transonic  wing/body 
calculations  in  which  the  full  potential  equation  is  solved  using  coordinates 
that,  in  general,  do  not  coincide  with  the  configuration  surface.  The 
wing/body  surface  appears  as  a  curved  boundary  within  a  rectangular  computing 
grid.  In  order  to  enforce  the  f low-tangency  condition  exactly  on  this  curved 
boundary,  an  imaging  scheme  is  used  which  involves  surface  control  points  and 
surface-adjacent  grid  points.  At  the  expense  of  some  computational  detail, 
the  approximations  inherent  in  a  small-disturbance  formulation  and  the 
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complicated  mappings  required  to  produce  a  surface-conforming  grid  are  both 
avoided.  Schemes  similar  to  the  one  outlined  have  been  applied  previously  to 
calculate  flows  about  simpler,  two-dimensional  planar  and  axisymmetric 
configurations. 

The  objective  of  this  contract  was  to  develop  a  computer  program,  based 
on  a  global  system  of  non-surface-fitted  coordinates,  to  calculate  transonic 
flowfields  about  wing/finite-fuselage  configurations  and  to  compare  calculated 
solutions  with  relevant  experimental  data.  It  was  determined  part  way  into 
the  contract  period  that  the  coordinate  formulation  originally  proposed  as  the 
basis  for  this  effort  was  inadequate  and  a  more  complicated  system,  involving 
a  higher  degree  of  nonorthogonality,  was  needed  to  achieve  acceptable 
solutions.  The  delay  associated  with  this  reformulation  as  well  as  some 
initial  difficulty  in  achieving  stable  operation  of  the  program  have  resulted 
in  a  geometry  capability  that  consists  of  a  general  wing  attached  to  fuselage 
of  semi-infinite  length. 

Sections  2.0  and  3.0  of  this  report  outline  the  formulation  of  the 
wing/body  problem.  Section  4.0  describes  the  numerical  solution  scheme.  A 
discussion  of  results  is  given  in  Section  5.0,  and  remarks  relating  to 
required  code  improvements  and  possible  future  directions  are  contained  in 
Section  6.0.  Appendix  A  is  a  summary  user's  guide  for  the  computer  program, 
and  Appendix  B  is  a  program  listing. 
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2.0  GEOMETRY  AND  COORDINATES 


Consider  the  configuration  shown  in  Figure  1.  A  wing  with  multiple 
sweep,  taper,  and  dihedral  is  attached  at  a  mid-height  position  to  a  blunt- 
nosed  cylinder  of  semi-infinite  length.  The  wing  leading  and  trailing  edges 
are  comprised  of  straight-line  segments;  each  edge  may  (but  need  not)  contain 
a  single  break  at  a  spanwise  location  that  differs  from  the  other.  The  wing 
varies  arbitrarily  along  its  span  in  section  profile  and  section  twist 
angle.  It  is  assumed,  however,  that  the  variation  of  wing  properties  is 
piecewise  linear  between  defining  sections.  The  fuselage  component  of  the 
configuration  has  an  asymmetric  side  profile  and  a  circular  crossplane  section 
of  varying  radius  along  its  length. 

For  convenience  in  formulating  the  numerical  scheme  to  calculate  the 
flowfield  about  this  class  of  configurations,  a  global  system  of  wing-adapted 
coordinates  is  introduced  followed  by  a  series  of  coordinate  stretchings. 

This  accomplishes  several  purposes.  The  wing/body  geometry  is  transformed  to 
a  standard  and  relatively  simple  form  to  ease  computational  bookkeeping. 

Also,  the  infinite  physical  domain  is  mapped  to  a  finite,  rectangular 


parallelopiped  along  whose  outer  boundaries  exact  far-field  boundary 
conditions  can  be  applied.  (However,  the  transformed  wing/body  surface 
remains  embedded  within  the  parallelopiped  as  an  irregular,  interior 
boundary.)  Finally,  grid-point  spacing  can  be  arranged  to  provide 
satisfactory  surface  resolution,  at  least  on  the  wing  at  the  present  stage  of 
development,  as  well  as  efficient  distribution  of  points  throughout  the 
computational  domain  (dense  near  the  configuration  surface  and  progressively 
more  sparse  as  distance  from  the  surface  increases). 

The  wing-adapted  coordinates  are  defined  by  the  following  relations: 


x  -  x,  (y) 


Y  =  X 

Y  b  * 


z  -  zm(x,y) 
T(y)  c (y ) 


there  x^e(y)  represents  the  streamwise  position  (sweep)  of  the  wing  leading 
edge,  c(y)  is  the  wing-section  chord,  b  is  the  wing  span,  T(y)  is  the  wing- 
section  thickness-to-chord  ratio  normalized  by  its  value  at  the  wing  root,  and 
zm(x,y)  is  a  nonplanar  mean  surface  determined  by  the  wing  geometry.  This 
mean  surface  is  specified  as  follows: 


V*.y)  *  *ie<y> 

[x  < 

xle(y)]’ 

(2a) 

=  *le(y) 

-  [x  -  xle(y)]  tan  at(y) 

^xie (y ) 4  x  * 

xte(y)1. 

(2b) 

"  Zle(y) 

-  c (y)  tan  <*t(y) 

[x  > 

xte(y)1. 

(2c) 

where  zie(y)  denotes  the  vertical  position  (dihedral)  of  the  wing  leading 
edge,  «t(y)  is  the  wing-section  twist  angle,  and  xtg(y)  gives  the  streamwise 
location  of  the  wing  trailing  edge.  Note  that  c(y)  =  x^fi(y)  -  xte(y). 

The  transformation  Equations  (1)  effectively  shear  out  wing  sweep,  taper 
in  both  chord  and  thickness,  dihedral,  and  twist.  In  the  (X,Y,Z)  coordinates, 
the  wing  appears  to  be  planar  with  a  rectangular  planform  and  uniform 
thickness  over  its  entire  span.  The  constant  in  the  x-trans format  ion  shifts 
the  X-origin  to  the  wing  mid-chord  locus;  the  leading  and  trailing  edges  are 
at  X  «  i  1/2,  respectively.  Lines  of  constant  X  coincide  with  constant- 
percentage-chord  stations  along  the  span.  The  Y  origin  is  at  the  wing/body 
centerline,  and  the  wingtip  is  at  Y  =  1/2.  The  Z-origin  is  on  the  mean 
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surface  zm(x,y).  Except  for  the  nonplanar  wing  mean  surface  zm(x,y),  the 
(X,Y,Z)  coordinates  are  similar  to  those  used  in  transonic  small-disturbance 
f ormulat ions . 

The  spanwise  region  containing  the  fuselage  and  the  region  outboard  of 
the  wingtip  are  treated  by  extending  the  wing  leading  and  trailing  edges  while 
holding  the  section  twist  angle  fixed  at  its  wing  root  and  wingtip  values, 
respectively.  The  fuselage  transforms  in  a  regular  but  not  well-defined 
manner  with  distortion  in  length,  profile,  and  crossplane  section.  In  dealing 
with  the  outboard  region,  the  normalizing  chord  length  in  Equation  (1)  is  held 
fixed  at  a  constant  value  c(y^)  beyond  an  arbitrarily  specified  spanwise 
position  y  =  y^  to  prevent  crossing  of  the  extended  wing  edges. 

In  order  to  produce  a  finite  computational  domain,  each  of  the  (X,Y,Z) 
coordinates  is  stretched  individually  as  described  below.  The  stretching 
formulas  are  adaptations  to  the  present  application  of  functions  used 
successfully  to  compute  the  transonic  flow  about  an  airfoil.  Figure  2  shows 
the  relationship  between  the  chordwise  and  spanwise  coordinates  and  also 
indicates  the  different  regions  of  the  physical,  intermediate,  and  stretched 
domains . 


•  X-coordinate  to  r-coordiuate 


Region  II:  X 

Region  I, III:  X 


f(at  +  a 2  <-2) 

*  X  +  A  tan  I” -p(r.  ±  f  ) 
o  1  2  o 


+  A.  tan 


[j(C  *  ro>3] 


(3a) 

(3b) 


In  Equation  (3b),  the  upper  set  of  signs  refers  to  region  I  and  the  lower  set 
to  region  III.  The  three  streamwise  regions  encompass  the  following  ranges: 


I:  —  <  X  v  -X  ,  -(1  +  t  )  <  f  <  -f  , 

o  o  o 


II:  -X  c  X  <  X  , 
o  o 


III:  X  <  X  <  °°  , 

o 


-r  <  f  <  r 

o  O 


<■  <  r  <  (1  +  r  ) 
o  o 


(4a) 

(4b) 

(4c) 
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Figure  2.  Coordinate  relationships. 


Regions  I  and  III  nominally  represent  the  regions  upstream  of  the  wing  leading 

edge  and  downstream  of  the  trailing  edge,  respectively,  while  Region  II  covers 

the  wing-section  chord.  This  stretching  is  symmetric  about  the  origin. 

Constants  a^  and  a2  are  determined  by  the  conditions  X  =  XQ  and 

dX/dF  =  irA,/2  at  £  =  F  .  The  transition  stations  ±X  between  the  cubic  and 
1  o  o 

tangent  stretching  functions  are  chosen  to  occur  a  small  distance  inside  the 
wing  leading  and  trailing  edges.  Parameter  F.  °  determines  how  much  of 
the  F -domain  is  confined  between  the  wing  edges.  For  an  evenly  spaced 
F-grid,  the  number  of  E -steps  in  Region  II  compared  with  the  total  number 
of  F -steps  will  be  in  the  ratio  C0/(l  +  Fq).  Constants  Aj  and  A£  mainly  allow 
control  over  grid-point  spacing  in  regions  I  and  III  but  also  provide  some 
control  over  the  uniformity  of  spacing  in  Region  II. 
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•  Y-coordinate  to  n-coordinate 


Region  IV:  Y  =  n  (b  ^  (5a) 

Region  V:  Y  -  Y^  +  tan  jj^(n  -  nQ)j  +  tan  ^y(n  -  n0^j  (5b) 

Spanwise  regions  cover  the  following  ranges: 

IV:  0  <  Y  <  Y  ,  0  <  n  <  n  ,  (6a) 

o  o 

V:  Y<Y<“,  n  <  n  <  (1  +  n  )  .  (6b) 

o  o  o 

This  spanwise  stretching  is  similar  to  the  one  used  in  the  streamwise 
direction  but  is  applied  only  to  the  half-space  because  only  unyawed  wing/body 
configurations  are  considered.  Region  IV  extends  over  the  wing/body  semispan, 
and  Region  V  encompasses  the  domain  outboard  of  the  wingtip.  Constants  bj  and 
b2  follow  from  the  requirements  that  Y  =  YQ  and  dY/dn  =  nBj/2  at  n  =  ho-  The 
transition  station  YQ  is  set  slightly  inside  the  wingtip,  and  parameter 
determines  the  fraction  of  the  n-domin  confined  to  the  semispan  region. 
Constants  Bj  and  B0  control  spacing  of  the  grid  outboard  of  the  wing  tip,  in 
region  V. 

•  Z-coordinate  to  £ -coordinate 

Z  «  Ci  tan  (ttc/2)  (7) 

Constant  Cj  controls  vertical  grid-point  spacing  near  the  wing  mean  surface. 
The  range  of  this  stretching  is: 

<  Z  <  »  ,  -1  <  C  <  1  .  (8) 

Figure  3  is  a  schematic  representation  of  the  nonorthogonal  (£,n,C) 
coordinate  system  which  indicates  its  relationship  to  the  wing  geometry. 

Figure  4  shows  a  schematic  representation  in  the  physical  (x,y,z)  domain 
of  a  grid-point  distribution  that  is  uniformly  spaced  in  each  of  the  (£,n,C) 


directions.  The  configuration  shown  involves  a  planar  wing  with  a  straight 
leading  edge.  Each  point  in  the  schematic  actually  represents  a  line  of  grid 
points.  Grid  taper  and  the  clustering  of  points  near  the  wing  mean  surface 
and  between  the  wing  edges  are  clearly  evident.  The  dashed  line  corresponds 
to  the  spanwise  station  y  =  discussed  previously.  The  planform  and 
crossplane  views  of  Figure  4  emphasize  the  wing  orientation  of  the  present 
coordinates  and  illustrate  their  lack  of  suitability  for  representing  the 
fuselage  region.  A  method  for  improving  resolution  about  the  fuselage  is 
described  in  Section  5.0.  The  wing-section  detail  emphasizes  the 
nonconformity  between  the  computation  coordinates  and  the  geometry  surface. 

For  later  use  in  writing  the  transformed  differential  equations  that 
govern  the  wing/body  flowfield,  the  stretching  functions  are  denoted 
symbolically  as 

4  =  i(X),  n  =  n(Y),  4  =  4  (Z ) .  (9) 

Stretching  derivatives  are  then  defined  bv  the  relations 

f(4)  =  d4/dX,  g(n)  =  dn/dY,  h  (4 )  *  dC/dZ.  (10) 

Evaluation  of  these  derivatives  follows  from  Equations  (3),  (5),  and  (7), 
respectively. 
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3.0  GOVERNING  EQUATIONS 


3 . I  Physical  Domain 

The  steady,  inviscid,  transonic  flow  past  a  wing/body  configuration  can 
be  considered  isentropic,  and  hence  irrotational,  if  only  weak  shock  waves  are 
present.  Such  a  flow  can  be  characterized  by  a  velocity  potential  $>(x,y,z) 
that  is  related  to  the  streamwise  x-component,  the  spanwise  y-component,  and 
the  vertical  z-component  of  velocity  by  the  relations 


u=it>  v  =  $  ,  w  =  .  (11) 

x  y  z 

where  coordinate-symbol  subscripts  denote  partial  differentiation.  In 
Equation  (11)  and  below,  it  is  assumed  that  all  velocities  are  normalized  by 
the  freestreara  speed  . 

Since  $  is  singular  at  infinity,  it  is  convenient  to  introduce  a 
perturbation  potential  $(x,y,z)  according  to  the  expression 

$  =  x  cos  a  +  z  sin  a  +  <p  ,  (12) 

where  a  is  the  angle  of  attack  of  the  incident  flow  measured  in  a  wing-section 
(x-z)  plane.  Then,  the  governing  equation  for  the  flowfield,  written  in 
cartesian  coordinates,  has  the  form 

(a2  -  u2H  +  (a2  -  v2)<(>  +  (a2  -  w2)<J>  -  2uv<j> 

xx  yy  zz  xy 

-  2vw<t>  -  2uw<f>  “0,  (13) 

yz  xz 

where 


u  “  cos  a  +  <p  , 
x 


(14a) 


V  =  ^  » 


(14b) 


w  *=  sin  a  + 


(14c) 
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The  local  speed  of  sound  a  is  determined  by  the  relation 


oo 


0  0  0  0 

where  q*  =  u  +  v  +  w  ,  M  is  the  freestream  Mach  number,  and  y  represents 
the  ratio  of  specific  heats  of  the  medium  (for  air,  y  =  1.4). 

The  pressure  coefficient  at  any  point  in  the  flowfield  is  given  by  the 
expression 


C 

P 


y  /  (y  —  l ) 


(16) 


The  boundary  condition  for  Equation  (13)  on  the  wing/body  requires  that 
the  flow  be  tangent  to  the  surface  and  can  be  written  as 


(cos  a  + 


V 


-  (sin  a  + 

x  y  y 


♦«) 


0; 


surface 


(17) 


derivatives  and  are  streamwise  and  spanwise  surface  slopes, 

X  y 

respectively.  In  addition,  the  Kutta  condition  requires  that  a  circulation 
T (yQ)  must  exist  at  each  spanwise  wing  station  yQ  which  is  of  such  magnitude 
that  the  flow  passes  smoothly  off  the  sharp  trailing  edge.  A  vortex  sheet 
extends  behind  the  wing  whose  strength  corresponds  to  the  spanwise  variation 
of  the  circulation.  Using  a  linearized  model  that  neglects  roll-up,  we  assume 
that  the  vortex  sheet  coincides  with  the  wing  mean  surface,  defined  by 
Equation  (2c),  between  the  trailing  edge  and  downstream  infinity.  There  is  a 
jump  r(yQ)  in  the  potential  function  across  the  sheet  which  is  constant  along 
lines  lying  in  the  sheet  that  are  parallel  to  the  freestream  direction.  The 
normal  component  of  velocity  and  the  pressure  must  be  continuous  through  the 
sheet,  however.  Thus,  on  the  vortex  sheet. 


i  (x  ,y  ,z  +)  - 
te  o  m 


♦  (».  ,y  .z,,  ) 

te  o  m 


F  (y  ),  41  continuous, 
o  z 


(18) 
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Far  from  the  wing/body,  the  flow  is  undisturbed  except  in  the  downstream 
Trefftz  plane  where  the  vortex  sheet  induces  a  two-dimensional  downwash 
flow.  Thus, 

4>  (*)  *  0  (19) 


everywhere  except  on  the  y-z  plane  at  downstream  infinity  where  the  potential 
function  satisfies  the  equation 


(a  2  -  v2)  <(>  -  2vw<$  +  (a  2  -  w2)  <j>  =0 

oo  Tyy  Tyz  oo  7.Z 


(20) 


subject  to  the  boundary  conditions 


<j>  </’  ~  (sin  a  +  <t>  )  = 

y  y  z 

L  J  J  surface 


=  o  , 


(21a) 


<j>(“,y,z^)  -  4>  (°°  »y » zm )  =  f(y)  on  the  Kutta  slit. 


(21b) 


<J>  ♦  0 


at 


/y2  + 


z2* 


(21c) 


The  Kutta  slit  corresponds  to  the  crossplane  profile  of  the  vortex  sheet  at 
infinity.  From  Equation  (15),  it  follows  that  a^  =  1/M^.  Boundary  condition 
(21a)  applies  on  a  cross-section  of  the  semi-inf initely  long  fuselage. 

For  the  unyawed  case,  only  half  of  the  configuration  need  be  considered, 
and  the  symmetry  condition  v  »  0  is  imposed  on  the  vertical  plane  y  ■  0,  which 
contains  the  fuselage  centerline. 

Equations  (13)— (21)  collectively  define  the  problem  for  the  wing/body 
flowfield.  The  two-dimensional  problem  for  the  downwash  field  in  the  Trefftz 
plane  is  embedded  within  the  overall  three-dimensional  problem. 
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3 . 2  Computational  Domain 

The  computational  problem  is  formulated  by  transforming  the  equations  of 
Section  3.1  into  (E,n,C)  coordinates.  Substituting  Equations  (1)  and  (9)  into 
Equations  (13)  -  (14)  and  using  the  symbolic  definitions  given  in  Equation 
(10)  yields  the  equation  for  the  perturbation  potential: 


*<fVf  +  B<8Vn  +  C<hV<  +  D*C  +  EV  +  FV  '  J' 


where 


+  <a2  -  v2)  C2«,n)  -  2uv  f«,  , 

c  (n) 


B  =  (a2  -  v2)  g(n)/b2  , 


(  2  2 
r  \,  2  2.  ~l.  .  2  2.  ,  (a  -  w  ) 

I  t  (n)c  (n) 


2uvG(E  ,n  )G(E  ,n [uG  (E  , n )  +  vG(E,n,o]  h(0, 


D  -  |  j^(a2  -  V2)  G(E,n)  -  -~yj  f(E)g(n), 


o  r  9  )  ~  —  vw  ”1 

-  I  (a  -  v  )  G(E,n,C)  -  uvG(E,n)  +  'T  (n  )c  (V)  j  *(n)h(° 


F  =  2  (a2  -  u2)  G^y~  +  (a2  -  v2)  G(C  ,n)G(€  ,n,C)  -  uv 


r«,n)G«.n)J  [von.-.)  +^7]  f«)h«) 
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Figure  S.  Schematic  diagram  of  computational  domain 
showing  boundary  conditions. 


L  = 


.Vy(C,n)g(n)/b, 


G(C,n)  //;<S,n)  +  G(5,n,c>  •'/  «,n)  -y^—y 


h(c). 


N  *  (cos  a)  -  sin  a. 


The  far-field  boundary  condition  remains  ^  =  0  on  all  edges  of  the 
computational  domain  corresponding  to  infinity  except  the  transformed  Trefftz 
plane  where  the  governing  equation  becomes 


BT(g*  )„  +  CT(h*  )  +  ET<f>  =  JT, 

n  n  C  t,  n? 


(25) 


with 


BT  *  (a2  -  v2)  g(n)/b2. 
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CT  = 


2  2  ~2  (a-  ‘  w2) 

-  v  )  G^l-K^C)  +— - — 


-  2vw 


T  (n)c  (n) 


G(l-K0,n,O 

T(n)c(n) 


h(C), 


ET=  ¥ 


(a2  -  v2)  G(l+T  ,n,C)  -  ,  t 

“  o’  ”  T(n)c(n) 


g(n)h(c). 


JT  =  | 2vwl (n )  - 


(a2  -  v2)  H( 1+C 


( .n  )  |  (w  -  sin  ot)  T(n)c(n), 


and  velocities  v  and  w  follow  from  Equations  (23b)  and  (23c) 
with  <(>£.  =  0.  Equation  (25)  is  subject  to  the  boundary  conditions 


(1+S  ,n) 

o 


+ 


G(l-K0,n,c)  .^(l+Co,n) 


1 

t  (n)c(n) 


h(c)* 


-  sin  a  j  “  0, 

)  surface 


(26a) 


4»  (Cte  ,n  ,0+)  -  <(>(Cte,  n,0  )  =  r  (n )  on  Kutta  slit. 


(26b) 


6  *  0  at  n  ■  l+n  or  C  =  ±1.  (26c) 

o 

On  the  wing/body  symmetry  plane,  the  condition  v  =  0  is  enforced  using 
Equation  (23b). 


3.3  Local  Streamline  Coordinates 

For  use  in  applying  upwind-biased  finite  differences  (Section  4.1),  the 
potential  equation  is  rewritten  in  coordinates  that  are  locally  aligned  with 
the  stream  direction,  which  is  denoted  by  S.  In  such  coordinates,  the 
principal  part  of  Equation  (13)  takes  the  form: 

(a2  -  q2)  4>ss  +  a2  (M  -  *sg)  -  0,  (27) 
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where 
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t (n )c (n  ) 


v  G(?,n)  +  ' 


c(n ) 


f(Oh(C), 


W  = 


2uvl (n )  +  v  H(F  ,n  ) 


u~  cosa  -  G(£,h)  (w  -  sin  a)  t(n)c(n) 


c(n) 


+  <2v 


uI(F ,n  )  +  wl  (n  ) 


+  V2  H(£,n,C)>  (W  -  sin  a)  r(n)c(n). 


The  Laplacian  becomes 


A*  =  Pi(fV?  +  V8Vc  +  VhVc  +  Vfn  +  Tl*cc  +  Vf?  +  Wr 


where 


1  +G2(C,n) 


c2(n) 


f(£). 


Qx  *  g(n)/b  , 


R1  " 


G2(£,n)  +  G2(4,n,c)  + 


t2(i)c2(ti) 


h(C), 


Sl  =  2G(C,n)f(Og(n)/b, 


T1  =  2G(C  ,n  ,C  )g(n)h(C  )/b. 


V2  "  2 


G(g,n) 

c(n) 


+  G(F,n)  G(C,n,q) 


f<C)h(c), 


(31) 
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W  =  H(f  ,n) 


cos  a  -  C(^,n)  (w  -  sin  a)  r(n)c(n)  c(n) 


+  H(F,n,c)  (w  -  sin  a)  x(n)c(n)» 


3 . 4  Transformation  Derivatives 

In  Sections  3.2  and  3.3,  various  quantities  appear  which  represent 
transformation  derivatives  between  the  (x,y,z)  and  (X,Y,Z)  coordinates.  These 
derivatives  follow  from  Equation  (1)  according  to  the  following  definitions: 


G(x,y ) 


H(x,y) 


I(y) 


G(x,y,z) 


3X 

3y 


32X 

3x3y 


3  Z 


3y 


H(x,y,z) 


I(x,y) 


G  (x ,  y ) 


3  2Z 
3x3y 


3  Z 
3x 


Ky) 


2 

3  Z 


3  y3  z 


C(F,,n) 


H(C  ,n) 


I(n) 


G(C,n,c)  =  Gj(n,4)  +  G2(F  ,n) 

H(C,n,c)  =  Hjfn.c)  +  H2(C,n) 


I(C,n) 


G(C  ,n ) 


I(n) 


(32) 


(33) 


(34) 


(35) 


(36) 


(37) 


(38) 


(39) 
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Evaluation  of  these  transformation  derivatives  in  the  computational  domain 
uses  the  one-to-one  correspondence  that  exists  between  grid  points  in  the 
(x,y,z),  (X,Y,Z),  and  (f,  ,0,0  coordinate  systems.  The  splitting  of  the 
functional  dependence  shown  in  Equations  (35)  and  (36)  actually  occurs  in  the 
(X,Y,Z)  system.  Transformation  derivatives  which  are  not  included  among 
Equations  (32)  -  (39)  either  are  zero  or  have  been  explicitly  evaluated  in  the 
equations  of  the  previous  sections. 
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4.0  NUMERICAL  SOLUTION  SCHEME 


4 . 1  Finite-Difference  Approximations 

The  type-dependent  finite-difference  concept  introduced  by  Murman  and 
CoLe^  is  the  basis  for  numerical  schemes  to  compute  steady-state  transonic 
flowfields.  Central  differences  are  used  to  approximate  the  potential 
equation  at  subsonic  points  of  the  solution  domain,  and  upwind-biased 
differences  are  used  at  supersonic  points.  Thus,  the  mathematical  character 
of  the  equation  is  properly  represented  as  it  changes  type  from  elliptic  to 
hyperbolic.  The  original  application  involved  the  transonic  small-disturbance 
potential  equation. 

In  order  to  apply  upwind  differences  to  the  full  potential  equation,  it 
is  necessary  to  take  into  account  the  misalignment  between  coordinate  lines 
and  the  velocity  vector  at  any  given  point  of  the  flowfield.  The  principal 
part  of  the  equation  is  recast  in  a  form  that  constitutes  an  effective 
rotation  to  the  local  stream  direction,  denoted  by  S:^ 

(1  -  M2 )<hss  +  (Ad  -  *ss)  =  0,  (40) 


where  M  =  q/a  is  the  local  Mach  number  and  the  quantities  <)>  and  Ad>  are 
defined  in  Section  3-3.  Then,  upwind  differences  are  used  to  approximate 
contributions  to  the  first  term  on  the  lefthand  side  of  Equation  (40),  and 
central  differences  are  applied  to  factors  associated  with  the  second  term. 


Jameson  has  shown  further  that  the  relaxation  procedure  for  the  finite- 

difference  counterpart  of  Equation  (40)  can  be  viewed  in  terms  of  a  damped, 

three-dimensional  wave  equation  involving  an  artificial  time.^  The  time 

dependence  arises  because  of  the  appearance  in  difference  formulas  at  each 

grid  point  of  both  a  new  solution  value  from  the  current  relaxation 

step,  and  an  old  value  <f> ;  .  ,  from  the  previous  relaxation  step.  Thus, 

i  » J  t  R 

timelike  terms  occur  implicitly  which  correspond  to  ‘t’gj.*  These  terms  play  a 
role  in  controlling  the  stability  of  the  relaxation  process  but  do  not  affect 
the  final  solution  since  they  vanish  as  the  process  converges.  Sometimes  the 
damping  inherent  in  the  finite-difference  analog  of  Equation  (40)  must  be 


augmented;  this  can  be  accomplished  by  adding  to  Equation  (40)  a  term  of  the 
form 


A?st  -  -  t 


U  A  t  /  U  V  ,  w  \ 

~  I  <f  +  —  ifi  +  a  ) 

q  Ax  l  q  xt  q  yt  q  zt  I 


=  -  £  At  [^^^t  +  Q2g(n)<Pnt  +  R2h(0<i)?t  ’ 


(41) 


with 


fjjlu 

2  2 

q  c(n)A( 


P„  = 


c(n) 


+  vG((,  ,n) 


,  _  f  (€)u  v 

'2  2  .  ...  b  ’ 

q  c(n)A^ 


R2  =~ 


iillu. 


q  c(n)A£, 


uG(n)  +  vG(S,n.O  + 


w 


T(n)c(n) 


In  Equation  (41),  the  value  of  parameter  £  is  arbitrary  and  can  be  adjusted  to 
control  the  amount  of  damping  augmentation. 


In  writing  difference  approximations,  central  differences  based  on  old 
values  are  used  to  calculate  first  derivatives  required  to  evaluate  the 
velocity  components  in  Equations  (23a)  -  (23c).  At  grid  points  where  the  flow 
is  subsonic,  Equation  (22)  is  used,  and  second  derivatives  are  represented  by 
central  differences.  Typical  formulas  are 
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(fv4)4  - 


_ 1 _ )(f  .  f  \  .  (n) 

2(A£)2  \\  i+1'j»k  i.j.k J  9i+l,j,l 


-  <£i+i,j,it +  2£i,i.k +  £1-i,J,k)  [s  *!"l!k +  (>  -  i)  »S"’,k 

(n+l)  j 


+  (f .  .  .  +  f  .  ,  .  .  .  , 

i»J'k  | 


(42) 


(&Vn 


I 


2(A  n)‘ 


(o  +  a  )  (A  ^  —  A  ^  ) 

lsi,j+l,k  gi,j,kJ  ^i,j+l,k  vi,j,k' 


(g  +  g  )  (<)>('n+1^  -  4>^n+l^  ) 

l8i,j,k  gi,j-l,kJ  V9i,j,k  Vi,j-l,k; 


(43) 


U 


1 

4A£,A  C, 


4> 


(n) 

i+1 , j ,  k+1 


(n)  (n+l) 

^i+l.j.k-l  +  *i-l,j,k-l 


.(n+l) 
vi-l , j ,k+l 


(44) 


In  Equation  (42),  the  relaxation  factor  w  has  been  incorporated  into  the 
difference  expression.  Also,  in  Equations  (42)  and  (43),  stretching  function 
values  required  at  half-step  points  have  been  evaluated  as  the  average  of  the 
values  at  grid  points  on  either  side;  the  averaging  is  a  second-order 
approximation  to  the  stretching  function  at  the  half-step  location.  At  grid 
points  where  evaluation  of  the  local  velocity  indicates  the  flow  to  be 
supersonic,  Equation  (40)  is  used,  and  second  derivatives  contributing  to 
H>ss  in  the  first  term  are  upwind  differenced.  For  example,  in  the  case  where 
velocity  components  u,  v,  and  w  are  all  positive. 


<fv?  - 


2  (AO 


(f  +  f  )  (2<|><'n+^  -  <p  ^  -  ^n+1*  ) 

1  i,j,k  ri-l,j,k'  V^i,j,k  ?i,j,k  M-l,J,k; 


i-l,j,k  i-2,j,k  i-l,j,k  i-2,j,k  J 


(45) 
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L  (n+l)  (n+1)  ,  (n+1)  .(n+1) 

Ai Ac  [yi,j,k  ~  ^i.j.k-l  +  ^i-1 , j ,k-l  “  *i-l,j,k 


(46) 


Derivatives  associated  with  the  $gt  terms  are  approximated  by  upwind 
differences;  a  representative  form  is 


At(f<?Ct)  "  2A£  (fi,j,k  +  f 


)  r($(n+1)  -  *(n) 

i-l,j,k  [_^i,j,k  9i,j,t 


-  (<p(n+1)  -<p(n)  )1 


(47) 


If  u,  v,  or  w  is  negative,  Equations  (45)  -  (47)  are  revised  to  reflect  the 
appropriate  upwind  direction.  All  second-derivative  contributions  to  the 
term  (A<p  -  <j>  )  in  Equation  (41)  are  central  differenced  in  a  manner  similar 

to  Equations  (43)  and  (44). 

4.2  Boundary  Conditions 

Since  the  present  formulation  produces  a  finite  computational  domain 
(Figure  5),  application  of  far-field  conditions  is  straight-forward  on  those 
faces  that  represent  infinity  and  on  which  $  ■  0.  The  boundary  condition  in 
the  Trefftz  plane  (downstream  infinity  face)  depends  on  the  circulation 
distribution  and  is  not  known  a  priori;  therefore,  it  must  be  calculated  by  an 
embedded  relaxation  procedure  which  solves  the  two-dimensional  problem  for  the 
downwash  field.  This  calculation  is  based  on  Equation  (25)  with  boundary 
conditions  (26a)  -  (26c);  central-difference  approximations  are  used  similar 
to  those  at  subsonic  grid  points.  The  Tref f tz-plane  calculation  coincides 
with  updates  of  the  circulation  distribution.  Finite  differences  in  the  C~ 
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direction  that  traverse  the  vortex  sheet  are  adjusted  to  incorporate  the 
potential  juny  across  it.  At  the  wing/body  symmetry  plane,  where  v  =  0, 
Equation  (23b)  is  used  to  compute  values  of  the  potential  function  at  image 
points  located  one  grid-step  beyond  the  computational  domain.  In  order  to 
compensate  for  coordinate  nonorthogonality  at  this  plane,  derivatives  are 
replaced  by  one-sided  differences  that  are  oriented  to  correspond  with  the 
skewed  intersection  of  spanwise  grid  lines.  The  present  symmetry-plane  scheme 
is  an  extension  of  one  used  by  Boppe^  and  has  proved  to  be  numerically  stable. 

At  the  wing/body  surface,  special  treatment  is  required  to  enforce  the 
flow  tangency  condition  given  by  Equation  (24).  A  scheme  is  used  that  is 
based  on  ideas  derived  from  previous  numerical  applications  of  non-surface- 
fitted  coordinates . ^  qhe  Neumann  form  in  Equation  (24)  is  replaced  by  an 

equivalent  Dirichlet  condition  that  is  applied  at  image  points  below  the 
wing/body  surface.  Referring  to  Figure  6,  point  S  is  a  typical  surface 
control  point  defined  by  the  intersection  of  a  vertical  grid  line.  Using  the 
corner-shaped  boundary-point  array  shown,  one-sided  differences  based  on  the 
Lagrange  interpolation  formula  for  three  unevenly  spaced  points  are 
substituted  into  Equation  (24)  to  obtain  a  Dirichlet  value  for  the  potential 
function  at  point  S.  Extrapolation  along  line  S-l-2  then  transfers  this  value 
to  the  uniformly  spaced  image  point  D1  located  below  the  surface  boundary. 

Thus , 

<t> (D 1 )  =->  j//x(S),  //y(S),  <^1...4>6,  AS,  An,  AC,  ,  (48) 

where //^(S),  ,7^(S)  are  known  surface  slopes,  <f>  ^ . .  •  <f>  6  are  potential  values  at 
indicated  points  of  the  boundary  array  associated  with  point  S,  (AE,An,AO  are 
grid  stepsizes  in  the  respective  coordinate  directions,  and  6  is  the  offset 
distance  between  point  S  and  the  grid  point  above  it.  Repeated  applications 
of  the  Lagrange  formula  to  the  sets  of  points  connected  by  arrows  in  Figure  6 
provide  the  inter-grid  values  $...$.  The  potential-function  value  at  image 
point  D1  fixes  the  boundary  condition  during  relaxation  of  the  solution  at 
exterior  points  on  the  vertical  line  Dl-2.  After  each  relaxation  sweep  of  the 
exterior  field,  the  image-point  potential  value  is  recomputed.  In  the  event 
that  upwind  differencing  at  a  surf ace-ad jacent  point  necessitates  a  second 
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S  Surface  point  at  which  Neuman  flow-tangency 
condition  is  enforced 

1*6  Exterior  points  used  to  calculate  equivalent 
Dirichlet  potential  at  surface  point  S 

D1,  D2  Interior  image  points  used  in  numerical 
solution  scheme 

5  Grid-point  offset  distance 
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Figure  6.  Grid-point  array  used  to  numerically  enforce  the 
boundary  condition  at  a  typical  wing/body  surface  point. 


image  point  D2,  the  associated  potential  value  is  obtained  by  linear 
extrapolation  in  the  vertical  direction. 

In  order  to  treat  each  control  point  consistently,  the  corresponding 
boundary-point  array  is  usuflly  directed  to  the  exterior  side  of  the  wing/body 
surface.  Thus,  a  reorientation  of  the  array  occurs  at  a  locus  of  zero 
streamwise  slope  on  a  convex  part  of  the  surface,  as  shown  in  Figure  7.  This 
arrangement  has  the  effect  of  eliminating  direct  coupling  between  neighboring 
image  points  in  the  present  methodology.  The  boundary -point  array  is  not 
reoriented  in  concave  surface  regions,  however. 

Additionally,  it  has  been  found  useful  to  permit  surface  image  points  to 
be  multi-valued  in  order  to  more  accurately  represent  the  surface  condition  at 
configuration  edges.  Figure  8  indicates  schematically  the  situation  in  a 
wing-section  plane.  Image  point  P  represents  the  leading-edge  condition 
during  relaxation  of  the  solution  on  the  line  ILE-1  upstream  of  the  section. 
Then,  for  computation  on  the  segment  of  line  ILE  below  the  section,  point  P 
represents  a  lower-surface  condition.  In  the  same  way,  point  P  can  also  be 
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assigned  an  upper-surface  boundary  value.  Image  point  Q  at  the  section 
trailing  edge  is  treated  similarly.  Other  such  multi-valued  image  points 
occur  at  the  fuselage  nose,  along  the  side-edge  locus  of  the  fuselage,  and  at 
the  wing  tip. 

4.3  Computation  Procedure 

The  system  of  algebraic  difference  equations  for  the  wing/body  problem  is 
solved  iteratively  by  systemically  sweeping  through  the  computational  space 
(Figure  5)  and  relaxing  the  solution  along  vertical  columns  of  grid  points. 

In  the  present  methodology,  the  domain  is  swept  by  crossplanes  beginning  at 
the  upstream-infinity  boundary  and  proceeding  to  the  downstream-infinity 
boundary.  Each  crossplane  is  swept  by  vertical  lines  from  the  configuration 
symmetry  plane  to  the  spanwise-inf inity  face  of  the  domain.  In  a  crossplane 
that  intersects  the  wing/body  surface,  all  column  segments  beneath  the 
configuration  are  relaxed  first,  followed  by  column  segments  above  the 
geometry,  and  then  by  columns  in  the  outboard  region  of  the  crossplane.  The 
relaxation  procedure  is  continued  until  either  a  prescribed  convergence 
criterion  is  satisfied  or  a  specified  number  of  domain  sweeps  have  been 
completed. 

During  a  particular  iterative  cycle,  the  first  step  is  to  fix  the 
surface-boundary  condition  by  calculating  image-point  potential  values 
associated  with  each  surface  control  point.  The  solution  is  then  relaxed 
throughout  the  domain  as  described  above.  At  designated  intervals,  the 
circulation  distribution  is  updated  by  applying  Equation  (26b)  at  the  wing 
trailing  edge.  Following  each  circulation  update,  the  Tref f tz-plane  boundary 
condition  is  recalculated  by  an  embedded  relaxation  of  Equations  (25)  -  (26) 
for  the  downwash  field.  The  cycle  is  then  repeated. 

During  the  sweep  process,  image-point  values  for  the  wing/body  surface 
condition  are  substituted  sequentially,  as  required,  into  the  solution  array 
that  stores  the  potential  function.  This  procedure  simplifies  program  logic 
by  effectively  eliminating  the  distinction  between  surface-adjacent  grid 
points  and  interior  field  points  in  computing  finite  differences.  On  vertical 
grid  lines  that  intersect  the  wing/body  configuration,  grid  benchmarks  for  the 
surface-adjacent  points  are  used  to  designate  the  length  of  column-segments  on 
which  the  solution  is  relaxed. 
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5.0  RESULTS 


Representative  calculations  have  been  made  for  a  configuration  comprised 
of  ONERA  Wing  M6  attached  at  mid-height  to  a  hemisphere-cylinder  fuselage. 
Details  of  the  wing  geometry  are  given  in  References  18  and  19.  The  wing  has 
a  planform  with  a  leading-edge  sweep  angle  of  30°,  a  taper  ratio  of  0.56,  and 
a  uniform  section  whose  thickness  ratio  is  0.098.  The  fuselage  radius  is 
taken  to  be  25%  of  the  exposed  wing  semispan.  A  representation  of  the 
configuration  planform  is  shown  in  Figures  9  and  10. 

Figures  9  and  10  also  show  the  calculated  pressure  distributions  at 
several  spanwise  stations  on  the  configuration  for  a  freestr earn  Mach  number  of 
0.84  and  angles  of  attack  of  0°  and  3.06°,  respectively.  The  sharp  peaks  in 
the  wing  leading-edge  region  are  the  consequence  of  grid  coarseness  (see 
below).  There  is  an  indication  of  a  double-shock  structure  in  the  midspan 


Figure  9.  Calculated  surface-pressure  distributions  on  a  configuration  composed  of  ONERA  Wing  M6 
attached  at  mid-height  to  a  hemisphere-cylinder  fuselage:  Mx  =  0.84,  or  =  0°.  Grid  dimension: 

49  x  19  x  25.  Every  third  spanwise  station  is  shown. 
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Figure  10.  Calculated  surface-pressure  distributions  on  a  configuration  composed  of  ONERA  Wing 
M6  attached  at  mid-height  to  a  hemisphere-cylinder  fuselage:  Mx  =  0.84,  a  -  3.06°.  Grid  dimension: 
49  x  19  x  25.  Every  third  spanwise  station  is  shown. 

region  of  the  wing,  which  coalesces  to  a  single  strong  shock  near  the  wing 
tip;  this  effect  is  particularly  evident  for  the  lifting  case  in  Figure  10. 

In  both  cases,  the  pressure  distribution  at  the  fuselage  centerline  is 
modified  by  the  presence  of  the  wing,  again  with  a  more  pronounced  effect  in 
the  lifting  case. 

In  order  to  assess  solution  accuracy,  the  calculations  are  compared  with 
available  experimental  data  for  ONERA  Wing  M6  obtained  in  a  wing-alone  test 
(References  18  and  19).  Since  spanwise  grid  stations  and  test  stations  on  the 
wing  do  not  coincide,  the  computed  results  are  interpolated  linearly  along 
constant-percentage-chord  lines  to  the  positions  of  the  data  measurements. 

The  comparisons  are  shown  in  Figures  11  and  12.  In  general,  agreement  Is 
reasonably  good  especially  near  the  wingtip  where  the  presence  of  the  fuselage 
in  the  calculations  has  the  least  effect.  Differences  near  the  wing  leading 
edge  and  the  poor  resolution  of  the  double-shock  structure  along  the  wing  can 
be  attributed  to  grid  coarseness  in  the  computed  results,  while  trailing-edge 
discrepancies  are  more  probably  the  consequence  of  viscous  effects  in  the 
data.  (Note  the  definite  shock-induced  separation  exhibited  by  the  test  data 
at  the  wingtip  station  in  Figure  12.)  Anomalous  behavior  such  as  the 
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The  calculations  for  these  cases  were  carried  out  on  a  49  x  19  x  25  grid 
in  the  chordwise,  spanwise,  and  vertical  directions,  respectively.  Thirty- 
three  vertical  grid  lines  intersect  each  wing  section  chord,  thus  providing 
sixty-six  (upper  plus  lower)  surface  control  points  on  each  section  profile. 
Control  points  are  spaced  at  approximately  every  3%  of  section  chord  with  an 
edge-offset  of  0.5 %  chord  at  the  leading-  and  trailing-edge  points.  The 
pressure  peaks  in  Figures  9  and  10  correspond  to  the  second  chordwise  control 
point  from  the  leading  edge  on  each  section.  In  the  spanwise  direction,  13 
grid  stations  occur  on  the  configuration  semispan. 

The  calculations  were  performed  on  a  Control  Data  CYBER  175  computer  with 
FTN(0PT=2)  compiler.  The  nonlifting  case  of  Figure  9  required  2.10  min  of  CPU 
execution  time,  and  the  lifting  case  of  Figure  10  required  5.23  min. 
Convergence  is  based  on  a  10_,!t  cut-off  limit  on  the  maximum  correction  to  the 
potential  function  over  the  solution  domain  between  the  final  two  iterative 
sweeps.  Decreasing  the  cut-off  limit  to  10“^  Was  found  to  approximately 
double  solution  times. 

Table  1  shows  the  accuracy  of  the  scheme  used  to  numerically  apply  the 
surface  boundary  condition.  At  two  semispan  positions  on  the  wing,  velocity 
slopes  computed  from  a  converged  solution  are  compared  with  prescribed 
streamwise  surface  slopes.  Agreement  between  corresponding  slope  values 
extends  to  at  least  the  third  decimal  place,  one  order-of-magni tude  greater 
than  the  cut-off  limit  on  the  potential  function. 

In  order  to  explore  the  stability  of  the  computer  program,  additional 
calculations  were  performed  for  a  configuration  with  the  planform  shown  in 
Figure  4  in  which  the  wing  has  a  uniform  NACA  0012  section.  Convergent 
operation  without  the  use  of  damping  augmentation  was  achieved  in  nonlifting 
cases  for  freestream  Mach  numbers  up  to  0.99  and  at  selected  (supercritical) 
Mach  numbers  for  angles  of  attack  up  to  6°.  These  results  are  not  included  in 
this  report  since  no  comparison  data  exist  for  the  particular  geometry. 


*'  COMPARISON  OF  SURFACE-SLOPE  VALUES  AT  TWO  SEMISPAN  STATIONS. 
DZDX  =  PRESCRIBED  SLOPE;  SLOPE  =  CALCULATED  VELOCITY  SLOPE. 
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6.0  CONCLUDING  REMARKS 


The  present  work  substantiates  the  use  of  non-surface-fitted  coordinates 
for  numerical  computation  of  transonic  wing/body  flowfields.  However,  a 
number  of  modifications  and  improvements  remain  to  be  incorporated  in  order  to 
make  the  present  computer  program  useful  for  engineering  applications. 

The  capability  of  the  program  to  perform  calculations  on  a  series  of 
progressively  finer  grids  needs  to  be  completed  in  order  to  improve  computing 
efficiency.  This  modification  involves  two  requirements.  Verification  of  the 
grid-halving  subroutine  in  the  program  must  be  completed,  and  the  array  that 
stores  the  potential  function  must  be  moved  to  disk  storage,  with  a  sequential 
transfer  into  central  memory  of  only  those  array  segments  required  at  each 
stage  of  the  computation  process. 

The  chordwise  stretching  function  used  between  the  leading  and  trailing 
edges  of  the  wing  is  symmetric  about  the  mid-chord  locus.  In  order  to  improve 
resolution  of  the  blunt  leading-edge  region  within  a  fixed  grid  dimension,  the 
introduction  of  an  asymmetric  stretching  that  clusters  more  points  near  the 
leading  edge  than  the  trailing  edge  would  be  helpful.  This  modification  would 
provide  a  more  efficient  means  of  improving  solution  accuracy  in  the  leading- 
edge  region  than  a  simple  increase  in  the  number  of  chordwise  grid  points. 

Incorporation  into  the  program  of  a  finite-fuselage  capability  would 
provide  a  better  geometry  model  for  engineering  applications.  An  attempt  to  do 
this  during  the  present  contract  was  unsuccessful  when  the  resulting  version 
of  the  code  proved  to  be  nonconvergent .  The  principal  difficulty  seemed  to 
involve  the  question  of  how  to  properly  treat  the  vortex  sheet  (circulation 
distribution)  in  the  region  behind  the  fuselage. 

Finally,  there  exists  the  problem  of  improving  flowfield  and  geometry 
resolution  around  the  fuselage.  A  method  for  accomplishing  this,  depicted  in 
Figure  13,  involves  a  two-coordinate  formulation  in  which  the  present 
coordinate  arrangement  is  retained  about  the  wing  and  better-suited 
coordinates  —  perhaps,  a  simple  stretched,  cartesian  system  —  are  introduced 
in  a  vertical  slab  whose  width  coincides  with  that  of  the  fuselage.  This 
scheme  would  require  interpolation  within  an  overlap  region  (possibly,  only  a 
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single  plane)  common  to  both  coordinate  systems.  In  order  to  conveniently 
accommodate  the  two-coordinate  formulation,  it  may  be  preferable  to  revise  the 
computation  strategy  such  that  the  domain  is  swept  by  wing-section  planes 
beginning  at  spanwise  infinity  and  moving  to  the  wing/body  centerline.  The 
vertical-line  sweep  in  each  section  plane  would  begin  at  the  upstream  boundary 
and  proceed  to  the  downstream  boundary. 

The  baseline  computer  program  described  in  Appendix  A,  when  modified  as 
discussed  above,  is  expected  to  provide  a  framework  for  treatment  of  complex 
wing/body  configurations.  Further  extension  to  include  add-on  components  such 
as  nacelles,  stores,  or  additional  lifting  surfaces  should  be  possible  also. 


GP  I  ?  05J 7 8 

Figure  13.  Two-grid  arrangement  for  wing/body  flow 
analysis. 
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APPENDIX  A.  USER'S  GUIDE  TO  COMPUTER  PROGRAM 


A  computer  program  based  on  the  formulation  presented  in  the  main  body  of 
this  report  is  listed  in  Appendix  B.  The  subroutine  structure  of  the  program 
is  shown  in  Figure  Al,  and  a  summary  of  subroutine  functions  is  given  in  Table 
Al.  The  code  has  been  run  on  a  Control  Data  CYBER  175  computer  with  an  FTN 
compiler.  As  dimensioned,  it  requires  260k  (octal)  storage  locations  to  load 
and  execute.  The  present  form  of  the  program  does  not  use  peripheral  storage 
devices. 

The  solutions  presented  in  this  report  have  been  computed  on  a  single 
grid.  Some  logic  is  contained  in  the  code  to  permit  a  calculation  to  be  made 
on  a  series  of  progressively  finer  grids,  but  it  has  not  been  fully 
implemented.  In  particular,  subroutine  HALFS  interpolates  the  potential  array 
P(I,J,K)  from  an  initial  grid  onto  one  that  is  half-spaced  in  each  coordinate 
direction  of  the  computational  domain.  The  interpolated  array  is  then  used  as 
the  starting  solution  for  continuation  of  the  relaxation  process  on  the  new 
grid.  Input  parameter  NHALF  specifies  the  number  of  grid-halving  cycles  to  be 
performed.  Full  implementation  of  this  capability  will  require  that  the  array 
P ( I , J , K)  be  transferred  to  disk  storage  and  that  provision  be  made  for  a 
sequential,  plane-by-plane  transfer  of  P(I,J,K)  values  into  central  memory  to 
coincide  with  each  relaxation  sweep  through  the  computational  domain. 

Figure  1  defines  the  class  of  wing/body  geometries  that  can  be 
represented  by  the  program.  Input  data  are  smoothed  to  ensure  that  the  wing 
leading  and  trailing  edges  are  piecewise  straight  lines.  One  break  (kink)  is 
permitted  in  each  edge  but  need  not  be  present.  If  both  leading-  and 
trailing-edge  breaks  occur,  they  may  be  at  different  spanwise  positions. 
Although  the  formulation  outlined  in  Section  2.0  places  no  restriction  on  wing 
attachment  to  the  fuselage,  code  logic  assumes  that  the  horizontal  mid-plane 
of  the  grid  (5  -  0)  coincides  with  the  maximum  width  position  on  the 
fuselage.  Combined  with  the  specification  of  circular  fuselage  cross- 
sections,  this  computational  arrangement  effectively  limits  application  to 
mid-wing  configurations  at  present.  This  restriction  can  be  removed  by 
incorporating  a  two-coordinate  formulation  as  discussed  in  Section  5.0. 
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Figure  Al.  Subroutine  structure  of  the  'vlng/body  program. 
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TABLE  Al.  SUBPROGRAM  LIST 


NAME 

FUNCTION 

MAIN 

Reads  input  data;  controls  overall  program  logic 

HALFS 

Interpolates  solution  arrays  onto  a  half -spaced  grid.  (Has  not  been  verified.) 

GRID 

Defines  computational  grid;  calculates  wing  configuration  data  at  spanwise  grid  stations 

GEOM 

Calculates  wing/body  coordinates  and  slopes  at  planform  grid  stations  in  the  physical  domain 

ASPUNE 

Parameterizes  input  geometry  data  in  terms  of  arc  length  along  the  curve  prior  to  spline  interpolation 

SSPLINE 

Interpolates  input  geometry  data 

PROF  L 

Calculates  geometry-defining  quantities  in  the  computational  domain 

COEFF 

Calculates  fixed  quantities  which  appear  in  the  finite-difference  equations 

(NIT 

Inttiali2es  solution  arrays 

SOLVE 

Executes  one  relaxation  sweep  of  the  computation  domair ;  updates  the  circulation  distribution  and 
solution  boundary  values 

FARBC 

Computes  the  Trefftz-plane  boundary  condition 

SURFBC 

Calculates  boundary  values  at  control  points  on  the  wing/body  surface 

SURFSC 

(Entry)  calculates  boundary  values  at  side-edge  points  of  the  configuration 

SYMBC 

(Entry)  Calculates  image  values  at  the  wmg/body  symmetry  plane 

TRICOE 

Calculates  coefficients  of  the  finite-difference  potential  equation  at  grid  points  along  a  specified  vertical 
grid-line  segment 

TRIT 

(Entry)  calculates  coefficients  of  the  finite-difference  downwash  equation  in  the  Trefftz  plane 

INVERT 

Solves  the  finite-difference  potential  equation  along  a  vertical  grid-line  segment  to  obtain  updated  solution  values 

MAXI 

Determines  the  maximum  potential-value  increment  along  a  specified  vertical  grid-lme  segment  between  the 
current  and  preceding  relaxation  sweeps 

ARRAYS 

Prints  out  solution  arrays  to  specified  iteration  intervals:  circulation  distribution,  potential  distribution, 
surface  boundary  values  and  associated  image-point  values,  symmetry-plane  image-point  values.  (Used  mainly 
for  diagnostic  purposes.) 

CPCOMP 

Calculates  and  writes  the  pressure-coefficient  distribution  on  the  wmg/body  surface 

GP03-1190-1 


Table  A2  summarizes  the  sequence  and  format  of  input  data  required  by  the 
program.  Data  categories  are  as  follows: 

•  case  title, 

•  computational  grid  parameters, 

•  code  execution  parameters, 

•  case  specification  (Mach  number,  angle  of  attack), 

•  fuselage  configuration  data,  and 

•  wing  configuration  data. 

Suggested  values  for  grid  and  execution  parameters  are  included  in  the 
table.  Also  given  are  parameter-value  limits  imposed  by  current  code 
dimensions . 
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TABLE  A2.  GLOSSARY  OF  INPUT  DATA 


CARD 

COLUMNS 

VARIABLE 

EXPLANATION 

1 

1-80 

TITLE 

Case  title  (written  to  output) 

2 

1-10 

NXI 

Number  of  chordwise  grid  steps  at  start  of  calculation.  Maximum:  48 

11-20 

A1 

X  —  £  stretching  constant  in  Eq.  (3b).  Suggested  value:  0.16 

21-30 

A2 

X  —  £  stretching  constant  in  Eq.  (3b).  Suggested  value:  2.75 

31-40 

XIO 

Stretching  transition  point  of  £  —  coordinate;  see  Figure  2.  Ratio  XIO/  (1  +  XIO) 
determines  fraction  of  chordwise  grid  steps  which  occur  on  each  wng-section  chord 

41  50 

XCAPO 

Stretching  transition  point  of  X-coordmate;  see  Figure  2.  Must  be  inside  wmg-section 
edge.  Suggested  value:  0.495 

3 

110 

Number  of  spanwise  grid  steps  at  start  of  calculation.  Maximum:  18 

11-20 

21  30 

Y  -  17  stretching  constant  in  Eq.  (5b).  Suggested  value:  0.16 

Y  -  T?  stretching  constant  in  Eq.  (5b).  Suggested  value:  2.75 

31  40 

ETAO 

Stretching  transition  point  of  X)  —  coordinate,  see  Figure  2.  Ratio  ETAO(1  +  ETAO) 
determines  fraction  of  spanwise  grid  steps  which  occur  between  the  wing/body 
symmetry  plane  and  the  wingtip 

41-60 

YCAPO 

Stretching  transition  point  of  Y -coordinate;  see  Figure  2.  Must  be  inside  wmgttp. 

Suggested  value:  0.49995 

4 

1  10 

N2ETA 

Number  of  grid  steps  in  vertical  direction  at  start  of  calculation.  Maximum:  24 

If  20 

C! 

Z  -  f  stretching  constant  m  Eq.  (7).  Suggested  value:  0.45 

5 

1  10 

ITERM 

Maximum  number  of  iterations  to  be  executed  on  the  initial  grid 

1 1  20 

NHALF 

Number  of  grid-halving  cycles.  Set  NHALF  0  for  single  grid  calculation.  (Note: 

Subroutine  HALFS  has  not  been  fully  verified  ! 

21  30 

NPRINT 

Iteration  frequency  for  execution  of  subroutine  ARRAYS,  which  prints  complete 
solution  arrays  for  diagnostic  purposes.  If  NPRINT  0,  only  the  circulation  distribution 
is  printed  upon  completion  of  the  relaxation  procedure  If  NPRINT  >  ITERM,  complete 
solution  arrays  are  printed  after  the  relaxation  process 

6 

1  10 

WE 

Relaxation  parameter  for  potential  a»  elliptic  field  points.  Suggested  value  1.70 

11  20 

WG 

Relaxation  parameter  for  circulation.  Suggested  value:  1.00 

21  30 

DPLIM 

Convergence  cut  off  limit  for  maximum  potential-value  increment  between  successive 
iterations.  Suggested  value  10  4 

31  40 

EPSI 

Damping  factor  in  potential  difference  equation  used  at  hyperbofic  f<e(d  points. 

Suggested  value  0.00,  increase  if  instability  occurs 

7 

1-10 

ZMACH 

Freestream  Mach  number 

11  20 

ALPHA 

Angle  of  attack  of  the  wing  reference  plane  (in  degrees) 

8 

1  10 

NF 

Number  of  fuselage  coordinate  sets  to  be  read  from  the  following  cards  (one  XF,  ZF, 

RF  set  per  card) 

1-NF 

1  10 

XF 

x  coordinate  along  fuselage  beginning  at  nose,  see  Figure  1 

11  20 

ZF 

/  coordinate  of  fuselage  side-profile  reference  line;  see  Figure  1 

21  30 

RF 

Fuselage  crossplane  radius,  see  Figure  1 

9 

1  10 

JSECT 

Number  of  wing-section  data  sets  to  be  read  subsequently.  Maximum  5 

11  20 

JBL 

Sequence  number  of  wing  data  set  at  break  in  leading  edge.  If  no  break  exists,  set  JBL  -  0. 

21  30 

JBT 

Sequence  number  of  wing  data  set  at  break  m  trailing  edge.  If  no  break  exists,  set  JBT  0. 

The  fust  card  uses  alphanumeric  format  20A4.  All  remaining  cards  use 
repeated  floating-point  format  8F  10  5  Conversion  of  data  to  integer 
mode  is  performed  within  the  program  as  required 
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TABLE  A2.  (CONTINUED)  GLOSSARY  OF  INPUT  DATA 


CARO 

COLUMNS 

VARIABLE 

10 

1-10 

YS 

XLES 

ZLES 

1 

CS 

41-50 

ATS 

51-60 

TS 

61-70 

FS 

11 

1-10 

NSU 

11-20 

NSL 

21-30 

KSYM 

1-NSU 

1-10 

XSU 

11-20 

ZSU 

1-NSL 

1-10 

XSL 

11  20 

ZSL 

12 

- 

- 

13 

- 

- 

14-15') 

16-17  V 

_ 

_ 

18-19  J 

- 

- 

EXPLANATION 


Spanwise  station  of  wing  section 
x-coordinate  of  wing-section  leading  edge 
2-coordinate  of  wing-section  leading  edge 
Wing-section  chord  length 

Wing-section  twist  angle  relative  to  wing  reference  plane  (in  degrees) 

Wing-section  thickness-to-chord  ratio 

Repeat  indicator.  If  FS  =  0,  wing-section  coordinate  data  from  the  previous  span  station 
is  used;  the  next  card  specifies  parameters  for  the  wing  section  at  the  next  span  station 
(go  to  Card  12).  If  FS  =  1,  coordinates  for  a  new  section  profile  are  read  from  the 
following  data  cards. 

Number  of  wing-section  upper-surface  coordinate  sets  to  be  read  from  the  following 
cards  (one  XSU,  ZSU  set  per  card) 

Number  of  wing-section  lower-surface  coordinate  sets  to  be  read  from  the  following 
cards  (one  XSL,  ZSL  set  per  card) 

Wing-section  symmetry  indicator.  If  KSYM  -  0,  the  section  is  asymmetric;  both  upper- 
and  lower-surface  coordinates  must  be  given.  If  KSYM  =  1 ,  the  section  is  symmetric,  and 
only  the  upper-surface  coordinates  are  required. 

Upper-surface  coordinates  of  wing  section  from  leading  to  trailing  edge.  One  set 
per  card 

Lower-surface  coordinates  of  wing  section  from  leading  to  trailing  edge.  One  set  per  card. 
(Required  only  if  KSYM  *  0) 

Repeat  of  Data  Card  10  for  the  next  wing  section 

Repeat  of  Oata  Card  f  t  and  data  sets  f  -NSU  and  1-NSL  for  the  wing  section  defined  by 
Data  Card  12.  (Required  only  if  FS  =  1  on  Card  12) 

Up  to  three  additional  wing-section  definition  cards  and  coordinate  data  sets.  Total  number 
of  data  sets  must  correspond  to  JSECT  on  Data  Card  9. 


Fuselage  data  consist  of  sets  of  side-profile  reference-line  coordinates 
and  crossplane  radii  given  at  strearawise  stations  beginning  at  the  nose  and 
proceding  aft.  These  data  are  spline  interpolated  to  obtain  required  values 
at  grid-point  stations.  Along  a  constant-radius  segment  of  the  fuselage,  a 
number  of  values  must  be  given  to  ensure  accuracy  of  the  spline  fit,  and  the 
data  must  be  more  closely  spaced  near  the  ends  of  the  segment  than  in  its 
central  region.  A  final  table  value  at  streamwise  location  XF  >  25  is 
required  to  specify  the  semi-infinite  fuselage  length. 

Wing  data  are  read  section-by-section  beginning  at  the  wing  root  and 
moving  outboard  to  the  wingtip.  In  the  data  set  for  each  section,  the  first 
card  specifies  section  properties:  spanwise  position,  leading-edge 
coordinates,  chord  length,  twist  angle  relative  to  a  horizontal  reference 


line,  thickness-to-chord  ratio,  and  an  indicator  flag  to  designate  either  that 
new  profile  coordinates  are  to  be  read  or  that  the  previous  section  profile  is 
to  be  repeated.  The  next  card  specifies  the  numbers  of  coordinate  pairs  in 
the  upper-  and  lower-surface  data  tables  for  the  section  and  an  indicator  flag 
to  designate  whether  the  section  is  symmetric  or  asymmetric.  Next,  if  a  new 
section  profile  is  being  given,  the  upper-surface  coordinate  pairs  of  the 
profile  are  read  proceeding  from  leading  edge  to  trailing  edge,  followed  by 
the  lower-surface  coordinate  pairs  in  the  same  order.  If  the  profile  at  the 
new  section  is  identical  to  that  of  the  previous  section,  surface-coordinate 
data  are  read  internally  by  the  program  and  need  not  be  repeated  in  the  input 
data  list.  Also,  if  a  section  profile  is  symmetric,  only  the  upper-surface 
coordinates  are  required  in  the  input  data  set;  the  lower-surface  coordinates 
are  set  within  the  program. 

The  code  provides  the  following  tabular  output: 

•  input  data  listing, 

•  computational  coordinates,  stretching  derivatives  and  grid  benchmark 
values , 

•  wing  configuration  data  at  grid  stations, 

•  maps  of  surface-adjacent  grid  points  (upper  surface,  lower  surface, 
side  edge  in  planform  view), 

•  iteration  summary, 

•  final  circulation  distribution, 

•  detailed  solution  arrays  (if  requested),  and 

•  pressure-coefficient  and  surface-slope  distributions. 

The  input  listing,  coordinate  tables,  and  wing  data  provide  a  check  of  problem 
set-up.  The  grid-point  maps  define  the  wing/body  configuration  in  the 
computational  domain.  The  iteration  summary  printed  after  each  complete 
relaxation  cycle  lists  the  value  and  location  of  the  maximum  potential 
increment  between  the  current  and  previous  iterations,  the  number  ot 
supersonic  points  detected,  the  current  wing-root  value  of  circulation,  and 
information  about  the  embedded  iteration  process  required  to  update  the 
Tref f tz-plane  boundary  condition.  The  solution  arrays,  when  requested  via 
input  parameter  NPRINT,  include  the  complete  potential  array  as  well  as 
control-point  and  image-point  potential  values  which  arise  in  the  numerical 
application  of  the  surface  boundary  condition.  These  array  print-outs  are 
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useful  mainly  for  code  diagnostic  purposes.  Finally,  in  addition  to  the 
pressure-coef f icient  distribution,  both  prescribed  surface  slopes  and  velocity 
slopes  calculated  from  the  converged  potential  field  are  printed  to  provide  an 
accuracy  check  on  the  surface  boundary  conditions. 
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APPENDIX  B.  LISTING  OF  COMPUTER  PROGRAM 


PROGRAM  TWH3< INPUT. OUTPUT, TAPE 5* INPUT  *TAPl6*0UTPUT) 

SOLVFS  THt  FULL  POTENTIAL  EQUATION  FOR  TRANSONIC  FLOW  PAST  A 

general-wing/semi-infinite-fusllaoe  CONF iguration. 


c 


s 


c 


c 


DIMENSION  TITLE (20) *  X  TES ( 5 ) 


COMMON/CONST/ 

COMMON/OATAX/ 


COMMON/SOLVO/ 


ALPHA«ZMACH»OPLIM*EPSI » *E * WG *  BSP AN *  A I c * CS A  * RX 1  ,RX2* 

sna*tdeta,tdxi ,tdzeta,deta»oxi ,ozeta* ile* im, I  NOSE t 
ITAILtlTF, JM, JRO°T«OT IP»KM,XW* J81 * JB2* J83.PLCB1 * 

A5  .ETA0.NET  A,  NXI.NZETA.XCAPO  »X  10*  YCAPO* 
NF » XF ( 150)  *ZF ( ISO >♦ RE ( 150) ♦  JSEC  T » JBL  »  JBT  *  NSU ( 5)  * 


NSL (5) *  YS ( 5 ) ♦ XLES  <  5 ) *ZL£S<5) *CS(5>,ATS(b>*TS<5), 


FS (5 ) *  XSU ( 150 
~  ‘  [49*19'  “ 

)♦  I 


Bbi  U9* 


(150*5) *XSL(150* 

\  lBb^AX?3^MAX?fi 


,5) *Z5U< 150*5) 
<*9) , 
amma  ( 


*ZSL (150.5) 

m\mh" 


JMAXI*jMAXT«KMAXI*t^MAXT*NSUP.P(49*l9*^5)  »PLE1  ( 19)  * 
PNEW (25  I *PN0SE.PSYM(49.25)  . PT 1 ( 1 9. 25) . PT2 ( 1 5  *  25) ♦ 
PT3( 19,25)  , PTE  1 (19) *PT£2<19) ,PwBL (49*19) ,PW9U(49«19) 


<TITLE(N) *N=1,2 0) 

FNXI ,A1 *  A2  *  X 1 0  * XC apo 
FNETA, A3* A4,ETA0, YCAPO 
K NZET  A »  A5 


MST  ART  =  0 
NCYCLE=0 
RE AO  (5,901) 

READ  (5.902) 

REAO  (5*902) 

REAO  (5*902) 

NXI=FNXI 

neta=fneta 
NZETAsFNZETa 

read  (5.902)  FITERM*FNHALF*FNPRINT 
ITERMsFITERM 
NHALF  sFNHALF 
NPR INT  =  FNPRl N  T 

REAO  (5*902)  *E  ,  WG . DPl I M , EPS  I 
*  ZMACH. ALPHA 

(TITLE (N) ,N=1 *20) .ZMACH, ALPHA 
I  TERM,  NHALF,  NPR  I  NT  *  **£  *  WG  *  DPL  1 M ,  EPS  I 
NXT.A1*A2,XIO*XCAHO*NETA.A3*A4,ETAO. 


READ  (5*902) 


WRITE 
WRITE 
WRI  TF 


(5*9031 
(6,90*) 
(fa, 905) 


Wr  AHA  _  Ki  7  C  T  A  .  AW 


READ  (5.902)  FNF 
NFsFNF 

IF  (NF.EQ.O)  GO  TO  150 
00  110  Ns  1  *  NF 

110  READ  (5,902)  XF (N) *ZF (N) ,RF (N) 
WRITE  ( fa  *  9 1 1 ) 

DO  120  Ns  1 , NF 

120  WRITE  ( fa  *  9 1 2 )  N* XF (N)  * ZF (N)  *RF <N) 
150  CONTINUE 


READ  (5*902)  FJSECT.FJBL.F JBT 

JSECT  =F  JSECT 

J9l=FJBL 

JBTsf JBT 

J  =  0 

211  J=JM 


i 

! 


I 


i 


i 
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IF ( J.GT .JSECT )  GO  TO  237 

READ  (5,902)  YSU)  *XLES<  J)  .ZLES(J)  ,CS(J>  .ATS!  J)  ,TS(J)  »FS(J) 

IF  <FS( J)  . EO . 0  .  )  GO  TO  230 

REAL)  (5,902)  FNSU , FNSL * F*S YM 

NSU  < J) =FNSU 

NSL ( J) =FNSL 

KSYM=FKSYM 

NU=NSU  < J) 

00  218  1=1, NU 

REAO  (5.902)  xsu ( I . J> . ZSU ( I » J) 

218  CONTINUE 
NL=NSL  <J)r 

IF  (KSYM.EO.l)  GO  TO  220 
00  219  1=1. NL 

REAO  (5.902)  XSL < I  ♦  J>  * ZSL < I ♦ J) 

219  CONTINUE 
GO  TO  211 

220  00  222  1=1. NL 
XSL ( I . J) =XSU ( I . J) 

ZSL ( I . J) =  -ZSU ( 1 . J) 

222  CONTINUE 
GO  TO  211 
230  NSU ( J ) =NU 

DO  23S  1=1. NU 
XSU(I.J)=XSU(I,J-1) 

ZSU(I.J)=ZSU(I.J-1) 

235  CONTINUE 
NSL ( J ) =NL 

00  236  1=1. NL 
XSL(I.J)=XSL(I.J-1) 

ZSL(I,J)=ZSL(I,J-1) 

236  CONTINUE 
GO  TO  211 

237  CONTINUE 

WRITE  (6,921)  JSECT, JbL.JfiT 
00  380  J=l, JSECT 

WRITE  (6,922)  J , YS < J > * XLES < J ) , ZLES ( J ) ,CS ( J > » ATS ( J) * TS ( J) . FS ( J) 
IF  ( FS ( J ) .NE.O.)  GO  TO  370 
WRITE  (6,923) 

GO  TO  380 
370  NU=NSU(J) 

WRITE  (6,924)  NU 

WRITE  (6,926)  (XSU( I , J) ,ZSU( I ,J) .1=1 ,NU) 

NL=NSL ( J) 

WHITE  (6,925)  NL 

WRITE  (6,926)  (XSL(I.J) .ZSL(I.J) »I=1»NL) 

380  CONTINUE 

00  438  J=l» JSECT 

438  XTES (J) =XLES ( J) *CS ( J) 

IF  (JBL.Nt.O)  GO  TO  440 
JSECTX=JSECT-1 

XSLOP= ( XL L5< JSECT) -XLES ( 1 ) ) / ( YS ( JSECT ) -YS  (11) 
ZSL0P=(ZLES(JSECT)-ZLES(1) ) /(YS (JSECT) -YS(l) ) 

00  439  J  =  2 , JSEC  T X 

XLES ( J) =XLES (1) ♦ ( YS< J) -YS ( 1 ) ) *XSLOP 

439  ZLES ( J ) =ZLES ( 1 ) ♦ ( YS ( J ) -YS ( 1 > ) *ZSLOP 
GO  TO  445 

440  JBLX=JBL-1 

X5L0P= (XLES ( JBL ) -XLES ( 1 > ) / < YS  < JBL) -YS ( 1 ) ) 
ZSL0P*(ZLES(JBL)-ZLES(1) ) / ( YS ( JBL ) -YS ( 1 ) ) 

00  441  J  =  2  ,  JBLX 

XLES ( J) =XLES (  ) ♦(YS(J)-YS(l) )*XSLOP 
44)  ZLES  ( J) =ZLES ( 1 ) ♦ < YS ( J) -YS ( 1 )  ) *ZSLOP 
JBLP=JBL*1 
JSECTX=JStCT-l 

XSLOP= (XLES (JSECT) -XLES (JBL) ) /( YS ( JSECT ) -YS ( JHL ) ) 

ZSLOP= (ZLtS (JSECT ) -ZLES (JRL) )/( YS ( JSECT )- YS < JBL ) > 

00  442  JsJBLP. JSECTX 

XLES(J)=XLES(JBL) ♦ (YS(J)-YS(JBD ) *X$LOP 
44?  ZLES ( J) =ZLES< JBL) ♦(YS(J)-YS( JBL) )*ZSLOP 
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c 


c 


446  IF  (JRT.NE.O)  GO  TO  44  1 

XSL OP=  UTES  ( JaEC  T  )  -  xTtS  (  1  )  >  /  ( YS  ( USED  T  )  -YS  (  1 )  ) 

00  446  U  =  A . JSEC  T  X 

446  XTES(J)=XTES( 1 ) ♦ (YS(J)-YS(l) ) *ASLOP 
GO  TO  460 

447  JBTXsjHT-1 

XSLOP= (xTtS  < JdT ) -XTFS(1))/(YS(0BT> -yS ( 1 ) ) 

DO  44?)  J-d  »  JH  f  X 

44fl  XTES( J) =Xl£S< i) ♦ ( YS < J) -YS ( 1 ) )*XSLOP 
JATP  =  JBT  ♦  1 

XSLOP=  UTES<  JSECT)  -  XTfc  S(JRT) ) / ( YS ( JSECT ) -YS ( JBT ) ) 
00  449  JaJriTP* JSECTX 

449  XTES(J) =XTES ( JBT ) ♦ ( YS ( J ) -YS  < JB T ) ) *X5L0P 
45  0  CONTINUE 

00  461  U=1 » JSECT 
4S1  CS(J)=XTES(J)-XLES(J) 

BSPAN=2*YS (JStCT  > 


510 


520 


550 

59o 

600 


901 

90? 

903 


904 


905 


WPITF  <  6 ♦ 93 1 ) 

I  T  E  W 1  =  1 
GO  To  520 
wPITe  (6*932) 

ITEH1=1TE9M 

CALL  HALES ( I TEHM.NPRINT ) 
ITE«2=I TEM1*I TEWM-1 
CALL  GW  10 
GEOM 
PHOEL 

IN^T (MSTART) 

(6.941 ) 


CALL 
CALL 
CALL 
CALL 
WPl  Tp 


00  590  ITEHsI  TEW1 , 1TFX2 
CALL  SOLVE ( ITEM) 

IE  (NPWINT.EO.O)  GO  TO  550 

IE  ( F  TE R/NPW I N  T  *NpR INT.NE.1TEW>  60  TO  550 

CALL  ARRAYS (NPRINT > 

CALL  CPCOMP 

IF  (OPMAX.LE.OPLIM)  GO  TO  6 00 
CONTINUE 

CAlL  ARRAYS (NPPINT ) 

CALL  CPCOMP 

IF  (NCYCLb.EQ.NHALE )  STOP 

MSTART=1 

GO  TO  510 


FORMAT 

FORMAT 

FORMAT 


(20A4) 

(HE  10.5) 
(3GH1THANS0NIC 


WING/BODY  PROGRAM  TrfB3 


G.  E.  CHMIELEWSKI 
MCDONNELL  DOUGLAS  Rl 
ST.  LOUIS.  MISSOURI 
JUNE  1980 


FORMAT 


FORMAT 


52X43HCOOEO  BY: 

85X43H 
8SX43H 

85X43H 
/1X20A4/// 

5X40HMACH  NUMBER 
5X40HANGLE  OF  ATTACK 
3X7H0EGREES////) 

( 2 1H-E XECUT I ON  PARAMETERS// 
5X40HMAXIMUM  ITERATIONS 
5X40HGRIO  HALVING  CYCLES 
5X40HSOLUT JON  ARRAY  PRINT  CYCLE 
1A21HRELAXAT ION  PARAMETERS// 
5X40HELL IPT IC  POINT 
5X40HC IHCULAT ION 
5X40HCONVER6ENCF.  LIMIT 
5X40HOAMPING  FACTOR 
(22H-CU0H0INATE  PARAMETERS// 
5X4QHSTRF AMW ISE  DIRECTION 
5*40H 
5A40H 
5A4  OH 
5A40H 


[SEARCH 

63166 

LABS./ 

/ 

/// 

ZMACH 

g 

F8. 3/ 

ALPHA 

S 

F8.3. 

ITERM 

s 

18/ 

nhalf 

s 

18/ 

NPRINT 

* 

18//// 

RE 

g 

F  8  •  3/ 

"G 

s 

F  8 • 3// 

'PLIM 

s 

E12.3/ 

EPS! 

s 

F8.3/> 

NXI 

g 

18/ 

A  1 

3 

E  8  •  3/ 

A2 

3 

F8  •  3/ 

A  I  0 

S 

F8.3/ 

XCAPO 

3 

F8.3// 
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nro  nnnnnn 


« 

SX40HSPANWISE  DIRECTION 

NET  A 

s 

IB/ 

• 

5X40H 

A3 

3 

F  8 . 3/ 

5X4  OH 

A4 

5 

E  8.3/ 

♦ 

5X4  OH 

ET  AO 

3 

F8 . 3/ 

• 

5X40H 

Y  CAPO 

5 

F8.3// 

« 

5X40H VFRT I C AL  DIRECTION 

NZETA 

X 

18/ 

• 

5X40H 

A5 

3 

F  0 . 3 ) 

911 


FORMAT 

FORMAT 

format 


« 

* 

922  FORMAT 

* 

* 

* 


921 

924 

925 

926 
931 


FORMAT 

format 
format 
forma  t 
format 


» 

93?  format 

« 

* 

FORMAT 


941 


» 


(  in-///  1X1  3HFUSELAGE  DA T  A// /9X IHN ,  8X2HXF  ♦  8X2HZF  ,  8X2HRE/ ) 
(I10.8F  10.4) 

tlH-///inH  W 1  N(j  DAT  A/ / /20H  NUMBER  OF  SECTIONS  .110/ 

22H  SECTION  AT  L.E.  RREAK.  10/ 

22H  SECTION  AT  T.E.  HRtAK.Ib) 

(1H-///13H  WlNO  SECT  ION. 1 10// 

1X5HYS  =  F9.6.5XSHXLES=  FS),6.5X5HZLES*  F9.6, 

SX5HC5  =  F9.6.5X5HATS  =  F9.6.5X5HTS  *  F9.6, 

SXSHFS  s  F9.6) 

( 4 JH-5URF ACC  OATA  IDENTICAL  TO  PREVIOUS  SECTION) 

< 1 H-/26H  UPPFR  SURFACE  (X-2  PAIRS) // 1XSHNSU  =  IB/) 
<1H-/20H  LOWER  SURFACE  (X-Z  PAIRS) //1X5HNSL  *  18/) 

<  UFU«6> 

(1H-///4BH  INITIAL  START  ******************************* 

****##*• •••»»••»•»»»•« ) 

(1H-///45H  CONTINUATION  ON  HALF-SPaCED  GRIO  *4*********. 

51H«*««*«**»*«»»**«*»*****»****»«**«»«***«****««*»»*»»* 

40H*** ********************* ****************  ) 

UH-///lX17HlTtRATION  SUMMARY// 

6X4HI IFH. 10X5HDPMAX.4X1MI . 4 X 1 H J , 4 X 1 HK . 6X4HNSUP . 
14X6HGAMMAR.5X5HITERT. 7X6  iOPMAXT . 3X2HJT  » 3X2HKT / ) 


END 


SUBROUTINE  HALFS < ITERM.NPRINT) 

)  I RCUL AT  ION  ARRAYS  ONTO  A  HALF-SPACED 


INTERPOLATES  POTENTIAL  AND  Cl 
COMPUTATIONAL  GRID. 


COMMON/CONST/  ALPHA.ZMACH.OPLIMtEPSI »WE»wG.BSPAN.Al2»CSA.flXl »RX2» 
SNA.TD£TA.TOXI.TDZ£lA.OETA«OXI«DZETA*lL£»IMf I NOSE  » 
ITAIL»1TE»UM»UR00T.UTIP.KM.M«(.  J81 .  JB2*  J03.PLCB1  » 

COMMON/OATAX/  . A5.ETA0.NET A. NXI ♦NZETA.XCAPO.XIO.YCAPO. 

NF.XFU50)  »ZF(  150)  »«F(  150)  .  JSECT  »  JBL  ♦  J8T .  NS  J  (5)  » 

NSL ( 5 ) *YS(5) , XLES 15) ♦ ZLES (5) »CS <5) ♦ ATS (5) .TS (5) t 

FS(5)  *XSU(150.5)  .XSLU50.5)  »ZSy(i50.5)  »Zf  - 

COMMON/SOLVO/  DL1 <A9. 19) ♦DL2(A9»19> .DSi (49) *DS2(49: 

DU2U9.19)  »DPMAX»UPMAXT»GAMMA<19)  ♦  IM)  _  _  . 

«:BsSl»SV«4iTl: 

PT3( 19.25)  .PTE1 (19)  .PTE2(19)  .PW0L<49.  19)  .PWBU<49.19> 
DELL (49. 19) .DELS (49) .0ELU149. 19) .DZDXL (49. 19)  . 

OZDXU (49.19) .DZDYL (49. 19) .OZDYU (49. 19) .HWBl<49.19> . 
HWBU(49.19) »IZL(19) . I ZU ( 19) .JB00(49) .K80DL <49. 19) » 
KB0DU<49. 19) »ZL (49. 19) »ZU (49. 19) 


3  1  .  I  3 131  t 

ZSL(150.5> 

wini*" 


COMMON/SURF  / 


NXl»2.*NXl 
NFTA«2.*NETA 
n2£TA»2,*NZETA 
I  TERM* I  I ERM/2 
NPRINT*2*NPRINT 

IMNs2*IM-1 

JMN=2*JM-I 

KMN*?*KM-1 
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c 


c 


c 


IMNX=IMN-1 
JMNX=JMN-1 
KMNX=KMN-1 
DO  200  I=lfIM 
JE  =  JPOD (  I  ) 

JEX=JE-1 
DO  no  J=JE*JX 
DO  110  K  =  i  *  KM 
KN=KMN- (2*K-1 ) ♦ 1 
1 1 0  H(I fJfKN) =  P(  I* J,K) 

DO  120  K  =  2  t  KMNX  <  2 

120  P(I»J*K)=0.5*(P(I«J»k-1)*P(I,j»k*1)> 
130  CONTINUE 

IF  (JE.EQ.l)  DO  TO  20u 
DO  1P0  J=1 * JEX 
KE=KBOOJ ( 1  .  J) -1 
P(If J*KE)=0U1 (I.J) 

oo  no  k=kw,kh 

XN*KMN- (2*K-1 ) ♦  1 
140  P< I f JtKN) =P ( I fJ,K) 

KWN1=2*KW 

DO  150  K=K*N1 fKMNX.2 
150  P(l*jn)=0.5#(P(I,J.K-l)*P(l,j,K»l)) 
KE  =  KRODL ( 1  * J) ♦  1 
P< I *J»KE>  =DL1 < I« J) 

00  160  K*l*Krf 
KN=KMN- ( 2#K-  1  ) ♦  1 
160  P  < I f JfKN) =P ( I  *  J  t  K ) 

KWNl=2*KW-2 
DO  170  K  =  2tMrf'Mlf2 

1  7 0  P(I.J»K>=0.5*<P(I.JfK-1)*P(I,j»K*1)) 
1SI0  CONTINUE 
200  CONTINUE 


KWN=?«*w-l 
DO  210  1  =  1 NOSt f  I M 
JEXsJBUD ( I ) -l 
210  P  (  I  fJEXfKfiN)  =051  (I) 

PUNOSEf  IfKWNJ=PNOSF 
DO  220  J=1.JTIP 
P(ILE.JfK"N)=PLE1 <J> 

220  P<  ITE  .  JfKf«N)  =  PTE1  ( J) 

00  260  K=  1  f  KMhl 
00  250  J=1fJM 
DO  230  1=1. IM 
IN=IMN- (2*1-1 ) *\ 

230  P ( IN, J.K) *P< I f J,K) 

DO  240  I=2,IMnX,2 

24  0  P(I,J,K)=0.5*(P(I-l»jFK)^P(m.J»K)> 
250  CONTINUE 
260  CONTINUE 


310 


320 

330 

340 


P < 2* I  NOSE- 1 . 1 ,KWN> =DS1 (I  NOSE) 

P  <2*ILE»JTIPfKwN) =0S1 ( ILE) 
P(2*ITEfJTIP.KWN) =0S1 <ITE) 

no  340  k=i.kmn 

DO  330  I  *  1 » I MN 
00  310  J=1.JM 
JN=JMN-(2*JM-1) *1 
P<I,JN.K)*P(IfJ,K> 

DO  320  J=2,JMNXf2 

P(IfJ.K)=0.5*(P(I,J-1fK).P(1,J»1»K)) 

CONTINUE 

CONTINUE 


DO  410  J=1.JM 
JN  =  JMN-<2*J-1) ♦  1 
410  GAMMA <JN) =GAMMa (J) 

UO  420  J  =  2 . JMn X  ,  2 

420  GAMMA (J) =0 .5* (GAMMA ( J-l ) .GAMMA (U* 1 ) ) 
JT 1PN  =  2* JT IP-  1 
GAMMA (JT IPN.  1 ) =0. 


5! 


JR00TN=2*JR00T-1 
GAMMA  < JROOTn-2)  =0. 
RETURN 
END 


subroutine  grio 

CALCULATES  GRID-POINT  COORDINATES.  STRETCHING  DERIVATIVES.  AND 
RING  CONFIGURATION  DATA. 


DIMENSION  XI (49) . YCAP ( 19) 
COMMON/CONST/ 


COMMON/COORD/  CHORD <1 9). DCDY <1 9)  »L)TDY  ( 19)  .DXLEDY  ( 

n  r*  c*  m  .  r*  t  r  r  #  «  a  \  .  ri  a  la  i  c  .  r.  «  yxi  .  n  a  Ur  .  c: 


.FLE.FN.Ff 
T  AU ( 1 9 )  .X( 


COMMON/DATAX/ 


(19) . UTDY ( 19) .DXLEDY (19) ,ETA(19) ,F(4?> 
19) .GAMLE.GAMN.GAMF.GAMTE.GAMRT fH(25). 
9) . XCAP (49) .XLE(19) »Y(19) ,ZETA(25) .OF. 

^  i  nil  c  r\\j  i  %  r\  \  _  A  i  rtUi  t  /  i  fit  AuT 


I  AU  \  17  /  «  A  K**'*  9  l  V  I  \  f  f  T  \  fvr  i 

ZC AP  ( 25 )  .ZLE  (19)  .OZLEDY  (19)  *  ALPHAT  ( 19)  .DADY  ( 19)  »§wT 
A1 .A?. A3. A4.A5,ETAO*NETA.NXl.NZtTA»XCAPO,XlO»YCAPO, 
NF.XF (150) .ZF (150) .RF( ISO) . JSECT . JBL» J6T.NSU (5) * 

NSL (5) . YS  <  S ) .XLES(5) .ZLES(5) ,CS(5) .ATS (5) .TS(5)  ♦ 
FS(5)  .XSU(1S0.5)  .XSL<1S0.5)  .ZSUU50.5)  .ZSL  (150.5) 


10?  format 
1 
2 

ioa  format 

294  FORMAT 
1 
2 

355  FORMAT 

357  FORMAT 
1 

411  FORMAT 
2 

595  FORMAT 

gss^ormat 

1 

657  FORMAT 

734  FORMAT 
1 
2 

7749  FORMAT 

801  format 

036  FORMAT 


3 

4 

5 

6 
7 

871  FORMAT 
1 
2 


( 1M-///1X27HSTREAM*ISE  GRID  COORDINATES/// 

1XBHOXI  =  F14.6// 

9X 1 H I » 1 1X4HXCAP. 13X2HXI ♦ 14X1HF/) 

(I10»3(1X.F14. 6) ) 

(1H0/1X8HFLE  =  F14.6/1X8HGAMLE  =  F14.6/ 

1 X8HXC APLE  *  F 1 4 , 6/ 1 X8HX 1LE  =  F 1 4 . 6// 1 X0HFTE  a  F14.6/ 
IX8HGAMTE  =  F14.6/1X0HXCAPTE  a  F 1 4 . 6/ 1 X8HX I TE  *  F14.6) 
(//1X8HFN  a  F 1 4 . 6/ 1 X0HGAMN  *  F 1 4 . 6/ \ XBHXC APN  a  F14.6 

/1X0HXIN  s  F 1 4 , 6 ) 

(//lX8t-INOSE  =  I14/1XRHILE  =  I14/1X8HITE  *  114/ 

1 X8HI T  AIL  =  I14/1X8HIM  =  114) 

( 1H-///1X25HSPANWISE  GRID  COORDINATES/// 

1 A8HDET A  a  F 1 4 .6// 

9X1hJ. 1 1X4HYCAP. 12X3HETA. 14X1HG/) 

(1H0/1X8HGWT  =  F14.6/1XBHGAMWT  *  F14.6/ 

1X8HYCAPWT  =  F14.6/1X8HETAWT  =  F14.6) 

(//1X8HGF  a  F14.6/IX8HGAMF  a  F 1 4 . 6/ 1 X0HYC APF  a  F14.6 

/1X8HETAF  a  F14.6) 

(//1X8HJR00T  a  114/1  X8H J  TIP  =  I14/1X8HJM  a  114) 

( 1H-///1X25HVERT ICAL  GRID  COORDINATES/// 

1 X8HDZET  A  =  F 1 4 , 6 // 

9X1HK. 1 1X4HZCAP. 1 1X4HZETA. 14X1HH/) 

(//1X8HKW  a  I14/1X8HKM  =  114) 

(It).10(lX.F10.6)  ) 

(1H-///1X35HRING  PLANFORM  AREA  S  *  F8.4/ 

1 X35H4 I NG  ASPECT  RATIO  AR  a  F8.4/// 

1 X35HST  AT  ION  AT  L • E  •  BREAK  YB1  a  F8.4. 

10X5HJB1  a  14. 1 0X7HE  T Ati 1  «  F8 . 4 . 1 0X7HPLCB 1  *  F8.4/ 
1X35HSTATION  AT  T.E.  BRtAK  Y02  a  F8.4. 

1UX5HJH2  a  14, 10X7HETA82  *  F8 . 4 . 1 0X7HPLCB2  *  T8.4/ 

1 X35HCONST ANT -CHORD  STATION  YB3  a  F8.4, 

10X5HJB3  a  1 4 , 1 0 X7HE  T AB3  =  F 8 . 4 » 1 0X7HPLCB3  =  F8.4) 
(1H-///1X41HWING  CONFIGURATION  DATA  (PHYSICAL  DOMAIN)// 
7X1HJ.8X3HXLE.4X7HDXLE/UY.6X5HCHORD.6X5HDC/DY.8X3HTAU. 
6X5HDT/DY,8X3hZLEi4X7HDZLE/DY.5X6HALPHAT.6X5HDA/0Y/) 
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on 


C 


xcap-stre  rcHiNG 

HP  I  =2 • *AT  AN ( 1 . ) 

X  I  M= 1 . ♦ X  I  0 
0X1=2. *XlM/NXl 

I  M  =  NX  I  ♦  1 
IOsNXI /2*1 
in=xio/oxi 
ILE=I0-I0 
ITE=IO*IO 
XI < I0)=0. 

II  =  10*1 

00  so  1=11. 1M 

So  XI  {  I  )  =XI  (  1-1 ) *0X1 
DO  60  1=2.10 
IX=I0-1*1 

60  XI  (  IX) =XI  ( IX* 1 ) -OX  I 

ALC  1=0.3*  <3.*XCAP0/XI0-HPI*A1 ) 

ALC2=0.S*  <HPI*A1-XCAPO/XIO) / (XIO*XIO) 

WRITE  (6.102)  OXI 
XCAP< 1 >=-l.£*99 
F ( 1 ) =0. 

ILX=ILE-1 
00  110  I*i.ILx 
IF  (I.EQ.l)  GO  TO  109 
PP  =  XI { I ) ♦ A 1 0 
TN2=TAN (HPI«PP) 

PP3=PP**3 
TN3=TAN (HPI*PP3) 

XCAP ( I ) =-XCAP0*Al*TN2*A2*TN3 

F ( I )  =  1./ (HP  I* (Al* ( l.*TN2*TN2> *3. *A2*PP*PP* < 1 . *TN3*TN3> ) ) 

109  WRITE  (6.108)  I  .XCAP  (  I  )  ,XI  (  I )  »F  (  1 ) 

110  CONTINUE 

DO  120  I=ILE.ITE 
PP2=X I ( I ) **2 

XCAP(I)=X1 <I)*(ALC1*ALC2*PP2) 

F  ( I) =1./ (ALC1*3.*ALC2*PP2) 

WRITE  (6.108)  I.XCAP(I) ,XI ( I) *F ( I) 

120  CONTINUE 

XCAP ( IM) b1.E*«9 
F  (  I  M )  s  0  « 

itep=ite*i 

00  130  I=ITEP»IM 
IF  (I.EU.1M)  GO  TO  129 
PP  =  XI  ( I ) -XIO 
TN2=TAN(HPI*PP) 

PP3=PP**3 

TN3  =  T  AN ( HP  I *PP  3 ) 

XCAP( I ) =XCAP0*A1*TN2*A2*TN3 

F(I)=1./(HPI*(A1*(1 ,*TN2*TN2) *3 . * A2*PP*PP* ( 1 . *TN3*Tn3) ) ) 

129  WRITE  (6.108)  I . XCAP (  I  ) , XI  (  I  >  .F  ( I  ) 

130  CONTINUE 

IF  (XCAPO.NE.O.5)  GO  TO  200 
X  ile  =  -x 10 
X  I  TE  =  X 1 0 
GO  TO  290 
200  XIl=XIO*OXI 
210  PP  =  XH-XIU 

TN2=TAN (HPI*PP) 

PP  3=PP  *  *  3 
TN3=TAN(HPI*PP3) 

FUN=A1*TN2*A2*TN1-0.5*XCAPO 

F  UNPW  =  HP I  * ( Al* ( 1 . *TN2»Tn2) *3.*A2*PP*PP* ( 1 .*Tn3*TN3)  ) 
OEL=FUN/FUNPR 
X I2=X I 1-OtL 

IF  (OEL.LT.l.fc-10)  GO  To  2S0 
X  1 1  =  X  12 
GO  TO  210 
250  X I  TE  =  X  1 2 
XlLEs-XI2 
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nno 


C 


290  GAMl_F  =  ABS  l  X ILE-X I  (ILF-1)  ) 

GAMTE  =  ABS (XI ( I TE ♦ 1 ) - X  I T E ) 

PP  =  XI TE-XIO 
TN2=TAN(hHI*PP) 

PP3=PP**3 

TN3  =  TA;vl  (HHI*PP3) 

XCAPTE=XCAP0*A1»Tn?*A2*TN3 

xcaple=-xcapte 

FTE  =  1 ./ (HP  I* (Al* ( 1 .♦TN2*TN2) ♦ 3. *A2*PP*PP* ( 1 .♦TN3#TN3) ) ) 
FLE=FT£ 

WHITE  (6.294)  FLF  .  GAMLF. .  XC  APLE  *  X I LE .  FTE  ♦  GAMTE .  XC  APTE ,  X I  TE 


OCS=CS (2) -CS ( 1 ) 
l)XS  =  XLES(2)  -XLES(  1  ) 

DYS=YS(2)-YS<  i) 

IF (DYS.EQ.O. >  STOP  »FRROR  YS ( 1 ) = fb ( 2 ) • 

XREF=XLES ( 1 ) -YS( 1 ) *UXS/DYS 
OXS=OXS<OCS 

CREFsXLES  (  1 )  ♦CSU  )-YS(  1 )  *DXS/DYS-XRtF 
XCAPN=-XPtF/C«EF-0.5 
00  300  1=1. ILE 
OIFF=XCAP ( I ) -ACAPN 
IF  (OIFF.GE.O.)  GO  TO  305 
300  CONTINUE 
305  I NOSF=  I 

XI1=XI (INOSE-I) 

310  PP  =  X 1 1 ♦ X I U 

TN2  =  T  AN ( HP  I *PP ) 

PP3=PP**3 

TN3  =  TAN  <HPI*PP3) 

FUN  =  A 1 *TN2» a2*TN3~XC  APN-XC  APO 

FUNPR=HPI* (A  I* ( 1 ,*TN2*TN2) ♦3.*A2*PP*PP* ( 1 .♦TN3*TN3) ) 
OEL=FUN/FUNPM 
X  1 2  =  X 1 1 -DEL 
0ELT=ABS(UEL) 

IF  (OELT.LT.l.E-lO)  GO  TO  350 
X 1 1  =  X  1 2 
GO  TO  310 
350  XIN=XI2 

GAMN  =  ABS (XI ( I  NOSE- 1 ) -XIN) 

pp=XIN*XIO 

TN2=TAN(HPI*PP) 

PPJ=PP«*3 

TN3=TAN«HPI*PP3) 

FN=1 ,/<HPI* (Al*< I .♦TN2*TN2) ♦3.*A2*PP*PP*( 1 .♦TN3*TN3) ) ) 
WHITE  (6.355)  FN.GAMN, XCAPN. XIN 
FL=XF(NF)-XF (1) 

XCAPTs(FL-XREF) /CREF-.5 
DO  360  1*1.  IM 

IF (XCAP(  I)  .GT.XCAPT)  GO  TO  370 
360  CONTINUE 
370  ITAIL=I-1 

IF(XF(NF) .GT.25.)  ITAIL=IM 

WRITE  (6.357)  INOSE . ILE ♦ I TE » I T AIL ♦ IM 


ycap-stretching 


ETAM=I .>£140 
ueta*etam/neta 
jm=neta> I 
jtip=i.*etao/ueta 

ETA  <  1  ) =0. 

DO  400  J=2.JM 
400  ETA (J) =ET A ( J- 1 ) ♦OETA 

ALC3  =  0.5*  <3.*YCAPO/ETaO-HPI*A3> 

ALC4s0.5* (HPI«A3-YCAP0/ETA0) / (E  T  A0#ET  AO ) 
WRITE  (6.411)  DETA 
DO  420  J=  i . JTIP 
PP2=ETA ( J) »*2 

YCAP( J) =ET  A ( J) * ( ALC3*ALC4*PP2> 

G ( J ) sl./(ALC3»3.*ALC**PP2) 

WRITE  (6.108)  J.YCAP(J) .ETAIJ) .G(J) 
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4  20 


445 

450 


500 

510 


550 

590 


C 


460 

600 

605 

610 


650 


C 

C 

C 


CONTINUE 
YC  AP (JM) =1 


G(  JM)  =0. 

JT  IPP= JT  I P* ♦  1 
00  450  j=jtipp,jm 
IF  (J.EO.JM)  GO  TO  445 
PP*tT A ( J) -ET&U 
TN2=TAN(HPI*PP) 

PP3=PP**3 

TN3=TAN(HPI*PP3) 

YCAP(J) =YCAP0*A3*TN2*A4*TN3 

5IJ) *1./ tHPI* (A3* ( 1 .♦1N2*TN2> ♦ 3 . * A4#PP#PP* < 1 . ♦ TN3* TN3)  )  ) 
WRITE  ( 6 » 106)  J*YCAP(J>  » E  T  A ( J) » G  ( J ) 

CONTINUE 


IF  (YCAPO.NE.0.5)  GO  TO  500 

etawt=etao 

GO  TO  590 

£TA 1 =£T  AO*D£  T  A 

pp=etai-etao 

TN2=TAN<HPI*PP) 

PP3=PP**3 

TN3=TAN(HP1*PP3) 

FUN*A3*TN2*A4*TN3-0.5*YCAPO 

FUNPP=HP  I  *  ( A  3*  (  l.-»TN2*TN?>  *3.  *A4*PP*PP*  ( 1 

OEL=FUN/FUNPR 

£T  A2=E  T  Al-UEL 

IF  (OEL.LT . l.t-10)  GO  TO  550 
ET  A 1 =£ T  A2 
GO  TO  510 
ETAWT=ETA2 

G4MwT  =  ABS(ETA(JTIPM)  -ETAWT) 

PP=ETA*T-ETAO 

TN2=TAN(HPI*PP) 

PP3=PP**3 


.♦TN3*TN3) ) 


TN3=TAN (HPI*PP3) 

YCAPW  T  =  YCAP0*A3*TN2* A4*TN3 

GWTsl ./  (HP  I*  (A3*  (  1.*Tin2»TN2)  ♦3.*A4*PP*PP*  ( l.*TN3*TN3)  )  ) 
WRITE  (6095)  GwT.GAMwT.  YC APWT .ETawT 


YFS=RF (1) 

00  460  N=2.NF 
IF (YFS.LT.Rf <N))  YFSsWF(N) 
CONTINUE 
ycapf=yfs/bspan 

00  600  J= 1  *  JT IP 

OIFF=YCAP' JJ-YCAPF 

IF  (OIFF.GE.O.)  GO  TO  605 

CONTINUE 

JROOT=J 

ETAl=ETA(JMOOT-l) 

etaisq=etai*etai 

FUN=FTA1*( ALCi* ALC4*FT AlSU) -YCAPF 
FUNPRsALCJ*3.*ALC4*ETA1S0 

oel=fun/funpr 
ETA2sET Al-OEL 

IF  (OEL.lT. I .£-10)  GO  TO  650 
ETA]=ETA2 
GO  TO  610 
ET  AF *E  T A2 

GAMF=ABS(ETA(JROOT)-FTAF> 

GF«1 ./ (ALC3*3.*ALC4«F  TAF*ETAF) 
WRITF  (6.655)  GF . GAMF * YCAPF » ETAE 
WHIFF  (6.657)  JROOT ♦ JT I P . JM 

zcap-stretching 

ZETAM=1 ,0 

0ZETA=2.*ZETAM/NZETA 

KM=N7ETA*i 

KW=NZE  T  A/i* 1 
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882«5 


ZETA (KW) =0. 

KK=KW* 1 

DO  710  K  a  K  K  *  K  M 

ZETA  (K) =  ZETA  (K-l ) ♦OZETA 

DO  720  K=2*KW 

KX=KW-K ♦ 1 

ZETA(KX)sZETA(KX*l)-DZETA 
WRITE  (6,734)  DZETA 
ZCAP ( 1 ) =- 1 .E*99 
H  (  1 )  =0  • 

ZCAP(KM>s-ZCAP(l) 

H (KM) =0. 

DO  7S0  K=1»KM 
ZET  =  ZETA  <K) 

IF  (K.EQ.l .OR.K.FQ.KM)  60  TO  749 
TN1=TAN(hPI*ZET) 

ZCAP(K) =Ab»TNl 

N(K)  =1  .  /  (Ab* HP I*(  1  •  ♦T»Ml*TfV»l )  ) 

WHITE  (6,108)  K,ZCAP(K) ,ZET,H(K) 

CONTINUE 

WHITE  ( 6  ,  /  749  >  KW , KM 

WING  CONFIGURATION  DATA 

JSL=J8L 

JST=JBT 

IF(JSL.EQ.O)  JSL=JS£CT 
IF(JST.EQ.O)  JSTajSECT 
XTIPsXLES(JSECT) 

CTIP=CS(JSECT) 

SWlNGsHSPAN* (XTIP*CTIP-XREF) 

SWING=SWING- (XLES ( JSL) -XREF) *YS ( JSL) 

SWING  =  SWING- ( XLES ( JSL) -XREF  *XTIP-XREF) * ( YS ( JSECT) -YS ( JSL) ) 
SWING=SWING- (ALES( JST) ♦CS( JST)-XHrF-CHEF)*YS(JST) 

SWlNGaSwING- ( X T I P*CT I P-XLES ( JST ) -CS (JST) ) *  < YS ( JSECT ) 4YS ( JST ) ) 
AR=8SPAN*dSP AN /SWING 

DO  800  Ja  1 » JM 
Y(J)aflSPAN*YCAP(J) 

CONTINUE 

CALCIaYS(JSECM-YS(JSECT-l> 

CALC2aXLES (JSECT) -XLES ( JSECT- 1 ) 

CALC3aCS (JSECT ) -CS (JSECT- 1 ) 

IF (CALC3.EU.0. )  GO  TO  803 

YSET=YS( JSECT )-CS< JSECT ) *CALCI/CALC3 

YB3= ( YSET4 YS ( JSECT ) ) /2  „ 

XLEBla XLES (JSECT) ♦ ( YB3-YS ( JSECT > > *C ALC2/CALC 1 

XTE83  =  XLES(  JSECT)  ♦CS(JSt'CT)  ♦  (  C  ALC3*CALC2  >  *  ( YB3-YS  ( JSEC  T  )  J/CAIC1 
DO  802  J= 1 , JM 
C ALCa YB3- Y ( J ) 

IF(CALC,LT,0.)GO  TO  805 

CONTINUE 

JB3= JM 

YB3= 1 • E ♦ 99 

PLCB3*0,0 

ETA83=ETA( JM) 

GO  TO  827 
JB3= J- 1 

YCAP83=YB3/BSPan 
ETAlaETA  (  JB3)  4l)E  T  A 
PP=ETAl-ETAO 
TN2=T  AN ( HP I *PP  > 

PP3=PP**3 

TN3=TAN(HPI*PP3) 

FUNaA3*TN2»A4»TN3-YCAPB34YCAP0  „  „  ,  _ 

FUNPRaHPI* (A3» ( 1 ,4TN?»TN2> 4 3 , * A4*PP*PP* ( 1 . 4 TN3*T NJ ) ) 

DEL  aFUN/FUNPH 
£TA2aETAl-UEL 

IF  < DEL ,Lr«l«E-10)  GO  TO  8826 

ETAI=ETA2 

GO  TO  8825 
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8826 


8833 


8834 


835 

8835 


8036 


8838 


E  T  AB3  =  £  T  4  2 

PLC8  3=  (E TAD 3-ETA  ( JB3)  )  /OET 4 
CONTINUE 


8839 


TO  8836 


00  828  J=1«JM 

IF  ( JBL.EQ.O)  GO  TO  830 

CALC*YS(JBL)-Y<J) 

ificalc.lt.o.)  GO  TO  B 3 L 

CONTINUE 

J8l=0 

Y81=YS<1> 

GO  TO  832 
J01=J-1 

ybI=ys ( jbu 

CONTINUE 

YCAPB1=Y81/BSPAN 
IF  (J81.EU.0)  JSTAPT=JR00T-1 
IF  (JB1.NE.0)  JSTA«T=JB1 
ETAI=ETA (USTAhT) 

etaisq=etai*etai 

FUN=FTA1*  ULC3* ALC4«ETAISQ)-YCAP8i 
FUNPRa ALC  3*  3 •*ALC4*ETA15Q 
DEL=FUN/FUNPR 
ETA2=ETA1-0EL 

IF  (PEL.LT.l.t-10)  GO  TO  8834 
E  T  A 1 =E  T  A2 
GO  TO  8833 
ETAB1  =  ET  A2 

PLCB1= (ETAB1-ETA (JSTART) )/DETA 

DO  835  J=l.JM 
IF(JflT.LQ.O)  GO  TO  8835 
CALCaYS(JBT>-Y (J) 

IF (CALC.LT.O. >  GO  TO  8836 

CONTINUE 

JB2  =  0 

YB2*YS ( 1 ) 

GO  TO  837 
JB?= J- 1 
YB2=YS(JHT ) 

CONTINUE 

YCAPB2=YH2/BSPAN 
IF  (JB2.EO.0)  JSTAHTajROOT-i 
IF  (JB2.NE.0)  JST  ART  =  JB2 
ETAlaFTA ( JSTART) 

etais5=etai*e  r  a  i 

FUN=FTA1* (ALC3*ALC4*E  TAISO) -YCAPB2 
FUNPRa ALC 3*3.* ALC4*ETA1S0 
OEl*FUN/FUNPR 
ETA2=t  T  A  1 -OEL 

IF  (DEL.LT.l.E-10)  GO  TO  8839 
ETAlaETA2 
GO  TO  8838 
ETABPaET  A2 

PLCB2= (ETAB2-ETA ( JSTART) ) /DET  A 
WRITE (6,836)  SW ING t AR * YB 1  * JB I ♦ tT AB I ♦PLCB 1 ♦ 
I  YB2»JB2»ETAB2»PLCB2»YB3»JB3»ETAB3«PLCB3 

SLOPF= (YS<2)-YS( l) )/(XLES( 2 )-XLEb(l>  ) 

DO  845  J= 1 « JB 1 

XLE  <J)»XLESU)  ♦  (Y(J)-YS(l)  ) /SLOPt 

oxledy (j)=i/ Slope 

CONTINUE 

CALC1=YS( JSECI ) -YS< JSECT-I > 

CALC2=XLES( JSECT > -XLES < JSFCT-1 ) 

SLOPE 1SCALC1/CALC2 
JB 1 Pa JB 1 ♦ i 
DO  847  J  =  JH 1 P  »  JM 

XLE (J)=XLES( JSECT) ♦ (Y(J)-YS< JSECT) ) /SLOPE  1 

DXLEDY(J)* 1 /SLOPE  1 

CONTINUE 
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c 


8 5?  DO  853  JS  =  2 « JSECT 

IF(YS(JS) .GE.Y (J) )  00  TO  8S7 

853  CONTINUE 

DYS=YS (JSECT) -YS< JSECT- 1 > 

OZLES*ZtES ( JSECT) -ZLES(J SECT-1) 
CALC=DZL£b/OYS 

ZLEB3=ZLES(JSECT) ♦<YR3-YS(JSECT) )*CALC 
IF(Y(J> .GT.YR3)  GO  TO  854 
ZEE  < J) =ZLEB3*  < Y ( J) -YB3) *CALC 
OZLEDY(J) =C  ALC 
GO  TO  858 

854  ZLE ( J ) *ZLER3 
DZLEDY ( J) =0 


GO  TO  858 

857  C4LC=(Z(,ES<  jS>-Zt_ES  (JS-1  ) ) / ( YS (JS> -YS /JS-I ) ) 
ZLE (J)~Z LtS ( JS> ♦7Y< J)-YS(JS) ) *CALC 
OZLEOY  < J) =CALC 


858  CONTINUE 


CALC= (CS (2) -CS ( 1 ) ) / <YS(2) -YS  t  1 ) ) 

00  888  JaliJM 

IF<YS<1)  .LE'.Y(J)  )  GO  TO  862 
CMORO <  J ) aCS  < 1 > ♦ <  Y  <  J > - YS ( 1 ) >  *C ALC 
OCOY ( J) =CALC 
GO  TO  868 

862  00  863  JS=2»JSECT 
1F(YS(JS)  .GE.Y  (J)  )  GO  TO  867 

863  CONTINUE 

DYS  =  YS < JSECT ) -YS< JSECT- 1  > 

OCS=CS ( JSECT ) -CS (JSECT- \\ 

CALCaOCS/UYS 

CHORORJsX IEB3-XLER3 

IF (Y ( J) ,GT .Y83)  GO  TO  864 

CHORD (J) =tH0RU«3* (Y(J)-Y83)*CALC 

OCOY(J)=CALC 

GO  TO  868 

86  a  CHORD  (J)  =CHOROH.I 
OCOY ( J) *0 
GO  TO  868 

867  CALCs (CS(JS)-CS(JS-n ) /(YS(JS)-YS(JS-l) ) 
CHORD (J) sCS  <  JS) ♦ ( Y ( J) -YS( JS>  >*CALC 

OCOY ( J»  =CALC 

868  CONTINUE 


TROOT  =  Tb (I) 

DO  B70  Jsl» JSECT 

870  rs(J) ~ts<j) /thoot 

CALC«(TS<2)-T5<1) ) / ( YS ( 2) -YS ( 1 )) 

DO  878  Jsl.jM 

I F ( Y S ( 1 ) .LE.Y(J) )  GO  TO  872 
fAU(J)=TSU)  MY<  J)-YSU>  )*CALC 
OTDY ( J) sCALC 
GO  TO  878 

87?  1)0  873  JS=*2»JSECT 

IF ( YS ( JS) .GE.Y (J) )  GO  TO  877 
873  CONTINUE 

OYS=YS (JSECT )-YS( JSECT- 1) 

OTS=TS(  JSECT)  -TS(  jSf-CT-1 ) 
CALC=OT$/OYS 

T AU8 3= TS (JSECT) ♦ (YR3-YS( JSECT) ) *CALC 
IF (Y (J) .GT.YB3)  GO  TO  874 
T  AU ( J) =T  AUB3* ( Y ( J) -YR3) *CALC 
DTDY (J) =CALC 
GO  TO  878 
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onnno n 


a  T  AU ( J )  =  T  AUHJ 
DTDY ( J) =0 
GO  TO  878 

7  CALC=< TS ( JS) -  I  5 ( JS-1 ) )/(rS(JS)-VS(JS-l)  ) 

TAU(J)=TS(JS) ♦ (Y ( J ) - YS ( JS ) ) *C  ALL 
0  TO Y ( J) =CALC 
3  CONTINUE: 

C ALC= <ATS<2)-ATS<l))/<YS(2>-YS<i>) 

DO  888  J=1*JM 

IF  < Y5( 1 ) .LE.Y ( J) )  GO  TO  882 
ALPHA  T  < J) =ATS( 1 > ♦ ( Y  < J) -YS( 1 ) )  *CALC 
OAOY ( J) =C«LC 
GO  TO  888 

?  DO  883  JS=2. JSECT 

IF  (YS(JS!  .GE. Y (J)  )  GO  TO  887 
B  CONTINUE: 

oys=ys  (  jsect  )  -ys  <  jsfc  r-n 

DATS=ATb (USECT ) -ATS (JSECT- 1 ) 

CALC=DATS/DYS 

ALPhTB3  =  ATS(JSECT)*(Y83-YS(JSECT))  <*CALC 

IF <Y (J) .GT.Y83)  GO  TO  884 

ALPHA  T  (  J)  =ALPr*TR3*  <Y(J)-Yh3)*CALC 

OAOY ( J) =CALC 

GO  TO  888 

►  ALPHAT  (  J.' =ALPHTH3 
OAOY(J)=0 
GO  TO  888 

p  CALC= (ATS (JS) -ATS (JS-1 ) ) / (YS ( JS) -YS ( JS-1 ) ) 

ALPHAT (J)=ATS(JS> ♦ (Y (J)-YS(JS) )*CALC 
0 AD Y ( J ) =C ALC 
I  CONTINUE 

WRITF  ( b  «  87  1 ) 

00  89S  J=1.JM 

.WHITE  ( 6 ♦ 8  0 1 )  J.XLE  < J)  .OXLFDY ( J)  .CHORD (J)  .OCDY(J)  , TAU(J> .OTDY(J) 
I  »  ZLF ( J )  .UZLEUY ( J,  , ALPHAT (J)  ,DAUY  l  J) 

00  900  I  —  1  »  I  M 
DO  900  J=i.JM 

X ( I »J) =<XCAP ( I ) *0.5) *CHORD( J) .ALE IU> 

CONTINUE 

RETURN 

END 


SUBROUTINE  GEOM 

COMPUTES  WING/BODY  COORDINATES  Z<X,  Y)  AND  SURFACE  SLOPES  DZ/DX  AND 
OZ/OY  AT  GRID-POINT  STATIONS  IN  THE  PHYSICAL  DOMAIN. 

01  MENS  I  ON  SF  ( 150)  .  XSET  <  150) ,S0(1 §0)  »  RF§  U  SO)  .  ZFS  ( 150)  ♦  XFSP  ,<  150 )  » 

*  RFSP(150) «  ZFSP ( 1 50 ) ,XFSPP<150> .RFSPP(ISO) .ZFSPP(lSQ) 

COMMON/CONST/  ALPHA.ZMACH.DPLIM.EPSI.wE.WGtBSPAN.AlZfCSA.RXl.RXZ. 

*  SNA.TDETA.TDXI.TOZETA.OETA.OXI.OZETA. ILE» IM.INOSE. 

*  ITAIL.ITE. JM. JR00T ♦UTIP.KM.KW.JB1 »JB2» JB3»PLCB1« 


COMMON/COORD/ 


CRORO( lbf?OCDr CJ9) »0T0Y(19> .OXLEDY < 19) ,£T A ( 1 9 ) ,F<4 9) 
»FLE.FN.FTE.G(19>  *GAMLE»GAMN»GAMF.GaMTE.GAMwT.H(2S) ♦ 
T  AU ( 1 9 ) ,  X  (  49 ♦  1  9) ♦ XC AP ( 49) » XlE ( 1 9) , Y ( 1 9) , ZET A ( 25) ,GF . 
ZCAP (25)  . ZLE (191  .OZlEDY (19) .ALPHAT  < 19) .DADY ( 19) ,§WT 

AJt  »  A? »  A3 .  AA *  Ab.E TAO.NETAjNX I  .nZET A. XCAPO. XIQ.  YC APB . 
NF.XFU50)  »ZF<150>  »RF<t50> .JSECT ♦  JBL  » JBT . NSU (5) . 

NSL ( 5 ) . YS ( 5) . XLES ( 5 ) ,ZLES<5| «CS ( 5 ) . ATS ( 5) . TS ( 5) . 
FST5) .XSU1150.5) »XSL( 150.5) . ZSU <150 . 5 ) . ZSL < 1 50 ♦ 5) 
PELL(49.191.DELS(4?T.DELU(49.19) ♦ OZ0XL < 49, 1 9)  , 


COmmon/dat ax/ 


common/suhf 


ofoXUU$!?9)  jD^DYL  (  4$??W!dz5|3  <4^J ?9UhS§L?U,  19) 
HWBU(49, 19) ,IZL U9) .  IZUU9) ,JB00(49) .K80DL (49. 19) . 

W  rt  Am  i  i  i  ft  .  ft  v  .  i  <  ifi  t  \  '...lift  .  ft  t 


"  '  v  ■  •  '  »  *  '  '  .  *  a.  w  '  •  ’  '  '  •  a.  v  *  a  f  #  y  vvf  vv 

KflOOU (49.19) , ZL < 4V, 1 9) ,ZU( 49, 19) 
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non  non 


i 

24-V  FOKMAT  (  1H-///1X29HPLANF0PM-A0JACENT  GRID  POINTS 

1  //3XlHl,«xTHJbOO(I)/) 

25 1  FORMAT  (14,11b) 


00  50  I  *  1 » I M 
JB00 ( I ) = 1 
00  50  J=1*JM 
ZU ( I  *  J ) =0 • 

ZL(I.J>=0. 

DZDXUU  ,  J)  =0. 

DZUXL ( I  * J) =0  . 

OZDYU(I,J)=0. 

50  OZOYL < I  * J) =0. 

FUSELAGE.  CO-ORDINATES 

00  56  I  =  INOSE  # I T A IL 
K= I - 1  NOSE ♦ 1 
DO  54  J  = 1  * JWOOT 
DO  5?  1 1  = INOSE  « I TA IL 
K  1  =  1 1  - I NO^E ♦ 1 
XSET(K1)=A<I1,J) 

52  CONTINUE 

N0=ITAIL-IN0S£*1 

IFCXSEr (NO) .Gr.XF(NF) )  NO=NO-l 

CALL  SSPLlNE (NF.XF,ZF,N0,XSET,ZFS»ZFSP,ZFSPP,4) 

CALL  ASPLINE  (imF  ,XF  , RF  , SF  , NO »  XSt- T  » SO , RFS »  XF SP  , RFSP » XFSPP » RFSPP * 2 ) 

IF  (K.GT.NO)  6=N0 

IF (Y(J) .GE.RFS <K) )  GO  TO  55 

CALC=SORT (RFS (K) *RFS (K) -Y ( J)*Y ( J) ) 

ZU(I » J) =ZFS (K ) ♦CALC 
ZL(I,J)=ZFS(K)-CALC 
DRFDXF*RFSP (K) /XFSP<K) 

OZFOXF=ZFSP(K) 

DZDXU  ( I » J)  =DZFDXF-»RFS(K)  *DRFDXF/CALC 
OZDXL (  I  ,  J) =DZFDXF-RFS (K) *ORFDXF/CALC 
DZDYU ( I , J>  =-Y ( J) /CALC 
DZOYL ( I  * J) =Y ( J) /CALC 

54  CONTINUE 

55  J0OD ( I ) =  J 

56  CONTINUE 

DO  124  I  =  ILE  *  I TE 
124  J80U(I)=JTIP*1 

WING  COORDINATES 

DO  150  US=1,USECT 
NU=NSU ( JS) 

DO  140  Ns i * NU 
XF  ( N )  ssXSU  ( N  « JS  ) 

RF (N) =ZSU(N, JS) 

140  CONTINUE  , 

DO  144  1=ILE»ITE 

XSET ( I ) = (X ( I , 1 ) -xlF ( 1 ) ) /CHORD ( 1 > 

144  CONTINUE 

NO«ITE-ILfcM 
00  146  1  =  ILL  *  I TE 
XSET ( I-ILt*l ) =XSFT(II 
146  CONTINUE 

CALL  ASPLINE (NSU< JS) » XF « RF , SF , NO , XSE T , SO » RFS » XFSP » RFSP , XFSPP » RFSPP 
1  »2> 

00  145  1*1*  NO 
ZSU(ILE*I-1.JS)sRFS(I) 

XSU ( TLE ♦ I - 1  * JS) =RFSP ( I ) /XFSP  <  I  > 
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c 


c 


c 


148  CONTINUE 
ISO  CONTINUE 

DO  154  I=ILE.1TE 
JS£CTx=JSECT-l 
no  154  jS=1»JSECTX 
00  153  J= JROO  T  « JT I P 

IF ( YS ( 1 )  .tO.Y ( J)  .OR. (YS( JS> .LT. Y ( J) .AND. Y  < J)  . LE . YS  <  JS* 1 ) ) ) 
1  GO  TO 
GO  TO  153 

IS?  CALC= (ZSU  < 1 ♦ JS) -ZSU < I . JS* 1 ) ) / ( YS( JS) -YS ( JS*1) ) 

ZU  (I  . J) =ZSU <  I . JS) *CALC*  <  Y  < J) -YS ( JS) > 

ZU  t  I  . J) =ZU l I  *  J ) *CHORf) ( J ) 

CALC=  (  XSU  1 1  ♦  JS)  -XSU  (  I  *  JS ♦ 1 ) ) / <  YS  ( JS )  -  YS  ( JS*  1 )  ) 

OZOXU  < I  * J) =XSU ( 1 » JS) ♦ ( Y ( J) -YS ( JS)) *CALC 

153  CONTINUE 

1 5 4  CONTINUE 


1SS 


156 


1ST 


DO  ISO  JS= 1  * JSECT 
NL=NSL ( JS) 

DO  1 55  N= 1 » NL 
XF (N>  =XSL IN. JS) 

RF (N) =ZSL (N, JS) 

CON  T INUE 

00  156  1=ILE»ITF 

XSET  U )  =  (X  < 1 . 1 ) -XLE  < 1 ) ) /CHOHD  ( 1 > 

CONTINUE 
N0=ITE-ILE*1 
00  157  1=1LE.ITE 
XSET (I-ILt*l)=XSET (I) 

CONTINUE 

CALL  ASPLINE  (nSL  ( JS)  .  XF  .  RF  .  SF  .  NO  ♦  XSET  » SO  •  RFS  ♦  XFSP » RFSP  ♦  XFSPR .  RFSPP 
1  .3) 


DO  1 55  1=1. NO 

ZSL ( ILfc*I-l • JS) =«FS< I ) 

XSL  (  I LE  ♦  I  “  1  .  J3)=P,"SP  (1  )  /  XFSP(  I  ) 

158  CONTINUE 

160  continue 

00  164  I  =  iL£  « I TE 
JSECT  X  =  JSEC  T-  1 
DO  164  JS= 1 .JSECT X 
00  163  JsJWOOT.JTIP 

IF(YS(l).tO«Y(J), OR •  ( Y  S ( JS ) • L  T  #  Y ( J ) . AND • Y ( J )  . LE • YS ( JS* 1 ) > ) 
1  GO  TO  1 6c! 

GO  TO  163 

16?  CALC= (ZSL U . JS) -ZSL < I . JS*1 > ) / ( YS ( JS) - YS ( JS* 1 ) ) 

ZL ( I . J) =ZSL ( I . JS) *CALC*  < Y( J) -YS ( JS)  ) 

ZL  < 1 . J) =ZL ( I . J) *CHO«U ( J) 

CALC= (XSL ( I . JS) -XSL ( I . JS* 1 ) ) / ( YS ( JS) -YS ( JS*1 ) ) 

DZOXL ( I . J) =X5L ( I . JS) ♦ ( Y ( J) -YS ( JS) ) *CALC 

163  CONTINUE 

164  CONTINUE 


DO  17?  I  =  I LE  »  I  T E 
DO  172  J=JHOQT . JTIP 
C4LC*-OXLtOY ( J) /CHORD (J) 

CALC=CAlC- IXCAP ( I ) ♦ .5) *DCDY ( J ) /CHOWO(J) 
CALC=CALC*ChOhD( J) 

CALC=CALC*0Z0XU< I , J) 

IF ( J.FU. JUl ,0H. J.F0.JB2.0R. J.Ew. JTIP)  GO  TO  170 
OZOYU ( I . J) =C ALC ♦ (ZU( I  * J)-ZU< I » J*l>  >/ ( Y ( J)-Y ( J*1 ) ) 
GO  TO  17? 

I  7 n  uZDYi)  (  I  .  J)  =CALC*  (ZU  (  I  » J)  -ZU  ( I .  J-l )  )  /  (  Y  ( J)  -Y  ( J- 1 )  > 
1 t ?  CONT [NUt 

00  18?  I  =  I LE  *  1  Tf 
00  18?  JSJROOT.JTIP 
CALC=-DXLEUY (j) /CHQWO( J) 

CALC=CALC- ( XCAP ( I ) *,5) *DC0Y (J) /CHORD (J) 
CALC=CALC*CHOHO ( J) 

CALC=CALC#OZOXl ( I . J) 

I F ( J.E0.JH1 .OR. J.F0.JH2.0R. J.EO. JT IP)  GO  TO  180 
OZDYL ( I . J» =C ALC ♦ (7L ( I . J) -ZL ( I . J*1 ) ) / (Y ( J)-Y ( j«l) ) 
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GO  TO  182 

1<*0  02UYH  I*  J>  =CALC*  (7L  (I»J)-ZL<  I*J-1M/(V  <J)-Y  (J-J)l 
18?  CONTINUE 

wRITF  (6.249) 

00  2S0  1  =  1.  IM 

250  WRITE  <b»25l)  I • JH00  ( I ) 

RETURN 

END 


SUBROUT InE  ASPLINE (NO. XO. YD . 5D »N0 » XO* SO . Y0 ♦ XP . YP . XPP » YPP . NFL > 

PARAMETERIZES  CXD.YD)  DATA  IN  TERMS  OK  ARC-LENGTH  SO  ALONG  CURVE  * 
THEN  SPLINE  FITS.  XP , YP , XPP , YPP  ARE  DERIVATIVES  W.R.T.  SO. 

NK'L.EQ.l  FOR  OX/OSs  0  (INFIMTE  OY/DX  AT  LEFT  END) 

NKL.E0.2  FOR  0Y/DS=*1  (POSITIVE  DY/DX  AT  LEFT  END) 

NFL. ED. .1  FOR  DY/DS=-1  (NEGATIVE  OY/OX  AT  LEFT  END) 


DIMENSION  XD (NO) ♦ YO (NO) . XQ (NO) *  YQ (NO) *SD ( ND) .SO (NO) 
DIMENSION  XP(NO) .YP(NO> ,XPP(NO) . YPP(ND) 


1 


EPSOsl.E-lG 

NlsNO-1 


SOU)aO 
Hi  =0 

0X1  =  XQ(2)-XDU) 
OY1=YO(2)-Y(j(i) 

Cl-SORT ( OX  1 **2  +  DY I **2) 
SD(2) =C1 

IF (NO ,E0 • 2 ) RE  T  URN 
00  L  1=2. N1 
0X1*X0  ( I )  -XDU-]  ) 
DY1=YD(I)-Y0(I-1) 

0X2  =  X0U*1)-X0(n 
DY2  =  YDU*l>~YU<I) 
0X*XD(I*1)-XD(I-1) 

0Y  =  Y0(I-»1)-YD(I-I> 

C2=SQRT  (DX2**?-.DY2**2) 
C=S0RT(UX**2*DY**2) 

A= (DYi«DX-OXl*OY) /? 
H=4*a/ (C*C1*C2) 

HAV= (HI.H) /2 

0S=C1* ( 1 ♦ (C1/2«HAV) *«2/6> 
S0<I)=S0<1-1MDS 
C  1  =C2 


Hl=H 

CONTINUE 

DS*C1*  < !♦ (Cl /2*H) **2/6 > 

SO (NO) =S0 (NO- I ) *0S 

CALL  SSPLINE (ND.SD.XO.NO.SO.XO.XP.XPP, 1 ) 
CALL  SSPLINE (nD.SD.YD.NO.SO.YO. YP.YPP.NFL) 
RETURN 
ENO 


SUBROUTINE  SSPLINE (NO .XO.YO.NO. X.Y.YP.YPP.NSW ITCH) 

S°L I NE  FITS  NO  DATA  POINTS  IXD.YU),  ASSUMES  CUBIC  RIGHT  END. 
LOCAL  ARRAYS  01Y.02Y.Q3Y  MUST  BE  APPRuPR 1 ATELY  DIMENSIONED. 

NSWITCH.EQ.1  FOR  CALL  FROM  ASPLINE  WITH  OX/OS*  0 

NSW I TCH.EO. 2  FOR  CALL  FROM  ASPLINE  WITH  DY/0S**1 

NSWITLH.E0.3  FOR  CALL  FROM  ASPLINE  WITH  DY/US=-1 

NSwITCH.EQ.4  FOR  DIRECT  CALL  AND  CUBIC  LEFT  ENO 
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D  IMENS  I  ON  X|>  ( i*L)  >  .  Yl)  (Mu)  . X (NO)  *  Y  (  'lO)  , YP (NO)  ♦ YPP (NO) 
,)  [  MENTION  D1  Y  U  SU)  t)>2Y  (  ISO)  »03Y  U  ?0) 

EPS  1 1 =- 1 . t- 1  0 
EPSIPs-tPSll 
NIMi=NL)-l 
UX=X0<2)-AD ( 1 ) 

If (Ox.EO.O. )  GO  T n  3S 
0F=<Y0<2>-Y U( l) )/0X 
0  1  Y  < 1 )=.b 
1,?  =  ? 

00  TO  (  l»2t3»*>  *  NSW  1  T  CH 
f)2Y ( 1 >  =3* (DF-O) /OX 
00  10  6 

02Y ( 1 ) =3* (OF-l ) /OX 
GO  TO  b 

0?Y  1 1 ) =3* IDF  *1) /Ox 
GO  TO  b 
12  =  3 

OXl  =XD  (  3  )  -XL )U) 

DF1= ( YO (3) -YO (?) ) /OXl 
C  =  (DX1**2-0X**2)  /OXl 
(3=  ( ox  ♦Ox  i )  *  (ox  *?*o x  l )  /nxi 
01Y (?) =c/ti 

02 Y  (?)  =t>*  (0F1-DF)  /B 
OX=DX 1 
IJF  =  DF  1 

00  T  I=I2fNlMl 
0X1=X0  U ♦ 1 ) -XU ( I ) 

IF  (UX1.EQ.0.)  GO  TO  36 
OF  1  =  (  Y  D  ( I  ♦  1 )  - 1 D  ( I )  )/Oxl 

B=2* (0X*DX1 ) 

F=6« (OF 1-OF ) 

DENOMsB-DX*DlY ( 1-1) 

02 Y  ( I )  =  (F-0X*02Y 11-1)  1 /OENOM 

01  Y  (  I ) =0X1/DEN0M 

0XS0X1 

OFaQFl 

CONTINUE 

0X1  =  XU  (NO-1)  -X0  (N(j-2) 

C  ALC  = ( OX ♦ OX  1 ) * (OX*OX1*2) /OXl 
CALC1= (OX  1 *#2-0 x*# 2) /UX1 
OENOMs  (CAUC-Oi  Y  (Nf)-2)  *CALC1  > 

02 Y (NO-1 1=  ( F-02 Y (NO-2) *CALC1 > /OENOM 


NlM2=N0-2 
00  8  I=1»NIM2 


K  =  NO— I  —  1 

IF  INSWITCH. EG. 4. ANO.K.EQ.l)  GO  TO  8 
02 Y  <K) =02 v <K> -01 Y (K) *02 Y  1 K ♦ 1 > 

CONTINUE 

02 Y  < NO )  *  (  (OK*OKl )  *f)?Y  (NO-1 )  -DX*D2Y  (ND-2)  )/OXl 
IF (NSW  I TCH.NE.4)  GO  TO  9 
OX=xn (2) -AO ( 1 ) 

OX  1 =XQ (3) -AD  <  2 ) 

02 Y ( 1 )  =  < (0X40X1 ) *02 Y (2) -0X*D2Y ( 3) ) /DX1 
K  =  NO 
00  11 

0X1=X0(K*1)-XD(M 

0F1=(Y0(K*1)-Y0<K) >/OXl 

0 1 Y  <  K ♦ l ) ~UFl*UXl/6* (02Y (K) *2*02Y <«.♦!)  | 

03Y  ( K  ♦  1 )  =  ( D  2  Y  ( K  ♦  1 )  -l)?Y  (K)  )/0Xl 
CONTINUE 

01  YU  I  =OF1-OXI/6*<2*D2Y(1>  *02 Y(2>  ) 

0  3  Y ( 1  I =03Y ( 2  > 

IF (NSwlTCH.EU. 1)  GO  TO  16 

INTEPPOLAJING  Y 


00  IS  U=1«N0 
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DO  12  1  =  1  *N[) 

OX  =  XD < I ) -X ( J) 

IF(DX.6E.fcPSIl.AND.0X.LE.EPSI2)  GO  TO  13 
IF <0X.6£.tPSI2)  00  TO  U 

12  CONTINUE 
60  TO  37 

13  Y(J)=YU(I) 

YP ( J) =0 1 Y  (  I ) 

YPP ( I ) =02 Y ( I ) 

60  TO  IS 

14  0X=X < J) -XU ( I ) 

Y  ( J)  =YD  ( I )  ♦DX*MD1  Y  ( I)  ♦OX/  2*  (Q2Y ( 1)  ♦0X/3*0JY  ( I )  )  ) 
YP ( J) =D1 Y (I ) ♦UX* (D2Y ( I ) ♦0X/2*03Y (in 
YPP ( J) =D2Y  < 1 ) ♦OX*D3Y  t  I) 

15  CONTINUE 
60  TO  23 

INTEWPOLAT IN6  X 

16  CONTINUE 

00  22  J=  1  ♦ NO 
00  17  1  =  1  *N0 

0Y=Y0(1T*Y(J»  , 

IF (0Y.6E.EPSI1.AN0.0Y.LE.EPSI2)  6U  TO  1R 
IF (UY.GE.tPSI2)  60  TO  1 9 

17  CONTINUE 
GO  TO  38 

18  Y<J)sYD(I) 

X< J)sXD(I) 

YP  < J) =01 Y ( I ) 

YPP ( J ) =D2  Y  ( I ) 

GO  TO  22 

19  0X=-DY/01Y ( I ) 

20  YO=YO( I ) ♦OX* (U1Y ( I ) ♦DX/2* (02Y (I) ♦DX/3*D3Y ( I) ) ) 

IF(DY7GE^tPSl 1 .AN0.0Y.LE.EPSI2>  GO  TO  21 
YOP  =  D 1 Y ( I ) *DX* (D2Y  < I ) ♦0X/2*U3Y ( I ) ) 

OELX=-OY/YOP 
OX=OX*OELX 
60  TO  20 

21  X(J)*XO(I)*OX 

YP ( J) =01Y ( I ) ♦OX* (D2Y ( I ) ♦DX/2*D3Y  (1) ) 

YPP ( J) =02 Y ( I ) ♦DX*D3Y ( I ) 

22  CONTINUE 


CONTINUE 
Wf TURN 

WRITE  (6.100) 

WRITF  (6.101)  X  D ( 1 )  »  X  0 ( 2 ) 

STOP 

WRITE  (6.100) 

WRITE  (6.102)  I.XD(I) *  XO ( I ♦ 1 ) 

STOP 

WRITE  (6.100) 

WRITE  (6.103)  J.X (J) .XO (NO) 

STOP 

WRITE  (6.100) 

WRITE  (6.104)  J, Y (J) .YO(ND) 

STOP 

FORMAT  (/SX.16HSU8R0UTINE  SSPLlNE/) 

FORMAT  (/SX.21HERR0R  IN  INPUT  XU <  1 ) =£ 1 2 . 4 ♦ 5X , 6HXD (2) =E 1 2. 4/ ) 
FORMAT  ( /SX , I6HERP0R  IN  INPUT  I  =  1 5 »5X . 6HX0 < I ) =E 1 2.4  * 5X . 

I  8MX0( t*I)*F12.4/)  ,  „  „ 

FORMAT  ( /SX . 23HX ( J >  IS  OUT  OF  RANGE  J= I 5» SHX ( J) =E 1 2 ,4 »5X » 
l  7MXD(ND)=E12.*/) 

FORMAT  ( /SX , 23HY ( J)  IS  OUT  OF  RANGE  U= 1 5  * SX .5HY < J > *t 1 2. 4 ♦ 5X . 
I  7HYD(ND) »E 12.4/) 


64 


nrn  r>nr>  o  nrnnnn 


r 

L 


SOBROUT I NE  PROFL 


CALCULATES  WING/HODY  AND  GRID  DATA  IN  THE  XJ-ETA-ZETA  COMPUTATION 
DOMAIN# 


DIMENSION  ZETAL  <49.  19)  .ZETAUU9. 19) 

COMMQN/CONST/  ALPHA « ZM ACH . OPL I M, EPS I » WE , WG ♦ BSP AN> A 1 2 ♦ CSA *  AX i . RX2 , 

►  SNA.TDETA.TOXI ,TDZ£ TA.DETA.OXI «DZET  A» ILE ♦  IM. I NOSE* 

>  ITAIL* 1TE. JM. JROOT. JTIP.KM.XW. JB1.UB2. JH3.PLC61. 

>  PLC82  *PLCB3 

COMMON/COOHO/  CHORO< 19)  .OCDY  (19)  »DTDY  (19)  ♦0XLE0YU9)  .ETA l 19)  *F  (49) 

>  «FLE*FN.FTE»(3(19)  »G4MLe.GAMN»GAMF.GAMTE»GAM*T»H<25)  . 

1  TAU(19) .X<49.19> *XCAP(49) »XLE ( 19) «Y(19) .ZETA<25) »GF. 

►  ZCAp ( 25) *ZL£ (19) .OZLEDY (19) .ALPHAT (19) f0A0Y(i9) .GWT 
COMMON/D A TAX/  A 1 .A2.A3.A4. A5.ETaQ.NET A .NX  I, NZETA.XCAPO.XiO.YCAPO. 

►  NF.XF (ISO) *ZF(150) *HF(150) . USECT ♦ JBL ♦ UBT »N$U < 5) ♦ 

k  N$L ( 5 ) »  YS ( 5) .XL£S(5) »ZLES<5) »CS (5) » ATS (S)  .TS<5). 

•  FS(5) *XSU ( 150.5) »XSL (150.5) «Z$U< 150*5) .ZSL(150.5) 

COMMON/SUHF  /  DELL  149.19) »D£LS(49> fDELU<49»19) *DZDXL<49»19) . 

-  DZ0XU<4  9, 19)  .DZDYLU9.  19)  . OZDYU ( 49.  19)  .HWBL  (49.19)  . 

►  HwBU  <  49 .19) «IZL(19) .IZU(19) .UB00U9) .KBODL ( 49. 19) ♦ 

>  KBOOU (49.19) .ZL(49.19) .ZU(49»19) 


304 

format 

1 

2 

306 

format 

309 

format 

1 

459 

FORMAT 

1 

461 

FORMAT 

601 

FORMAT 

1 

602 

FORMAT 

( 1H-///1X28HSUWFACE-ADJACENT  OHIO  POINTS/// 
‘  -  '  •'  "  DOWN.  J  ACROSS)// 


1A33HUPPER  SURFACE  (I 
3A1H1 .3X10HKRODU (I , J) /) 

(14.(2016)) 

(///1X33HLOWFR  SURFACE 
lUHKBOOL  < I ♦ J)/) 

(1H-///1X33H STATIONS  OF 
bHllU(J) .9X6HIZL ( j) /) 

UH-///1U0HNEGATIVF  UPPER-SURFACE  ghio-point  offset. 

24H  AT  PLANFOHM  STATION  1  =  I4.SH.  J  a  14/) 

( 1H-///1X40HNEGATIVE  LOWtW-SU«FACfc  GRID-POINT  OFFSET. 
24H  AT  PLANFORM  STATION  1  =  I4.5H,  J  *  14/) 


(I  DOWN,  J  ACROSS) //3XIHI.3X. 
ZERO  STREAMS ISE  SLOPE// 3X1  HO. 9X» 


100 


120 

C 


?DO 

C 


2  70 


HPI=?.*ATaN ( 1 . ) 
DO  100  1  =  1. IM 
DO  100  Jal.JM 
ZETAU ( 1 .0) 50. 
ZETAL ( I . J>  =0. 

DO  120  1  =  1.  IM 
DO  120  J=1»JM 
TC  =  T  AU  < J) *CHORD 
ZH=Zu ( 1 . J) /TC 
Z6  T  AU< I . J) *ATAN 
2  H  =  ZL ( 1  * J) /TC 
ZETAU  1  *  J  >  =  A  T  AN 


(0) 

(ZH/A5) 

(7H/AS) 


/HP  I 
/HP  I 


ah=as*hpi 
00  200  1=1 
DO  200  J=i 
ZH  =  ZE  T Au  ( I 
TN1=T AN (HP 
H*BU ( I . J ) s 
Zh=ZFTAL ( I 
TNI=TAN (HP 
HWBL ( I ♦ J) - 


.  I  M 
.JM 
*  J  I 
1  *ZH> 

1  ./  <  AH* 
.  J ) 
I«ZH) 

1 ♦/ ( AH* 


( 1 
( 1 


♦TNI* TNI ) } 
♦  )N1 *TN1 )  ) 


00  270  1  =  1. IM 
DO  270  J=I.JM 
KHOUU ( I .J) =0. 
KRODL (  I  . J) =0  . 
DU  30  0  1  =  1.  I M 
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on  300  J=1.JM 

Zh=Z E  TAU 1 1 » J) 

IF  (7H.EQ.0.)  GO  TC 
00  280  K=KW.«M 
IF  (ZH.GE.7tr4 (K ) ) 
KBOoU  < I »  J) =K 
GO  TO  28b 
CONTINUE 
ZH*ZFTAL( l*J> 

IF  (ZH.EO.O.)  GO  TO 
00  290  K=1»KW 
KK*KW-K*  I 

IF  (2H* LE.ZET A (KK) ) 
KROOl  <  I . J) =KK 
GO  TO  295 
CONTINUE 
CONTINUE 

continue 

WHITE  ( 6  *  304 ) 

00  JOS  I  =  A  •  1 M 
WHITE  ( b  * 306 )  I.<KE 
WHITE  ( b  *  309) 

DO  310  I  =  1 . 1 M 
WHITF  (6.306)  I. IKE 


I.  (KRODU  ( I  «U>  tJ~l  »UM> 
I  .  (KBOOL  ( I  *  J)  • J*1 » JM) 


00  400  J= 1 » JM 
IZU(vi)  =0. 

1ZL  <U) =0. 

nfx=nf- 1 

00  365  N=1.NFX 

IF  (RF(N*1) .LE.HF(N) )  GO  TO  370 

CONTINUE 

XF 1 =  XF ( N) 

00  375  l  =  INOSt . I  TAIL 

IF  (X (I .1 ) .GE.XF1)  GO  TO  380 

CONTINUE 

I F  S=  I 

JR00TX=JRU0T~1 

00  410  J=i.J«UOTX 

IZU(J) =IFS 

IZL  ( J) =lFb 

00  450  J=JROOT.JTIP 

00  420  I=1LE.1TE 

IF  (OZOXUU.J)  .GT.O.)  GO  TO  420 

IZU ( J) *1 

GO  TO  425 

CONTINUE 

00  430  I  =  1LE ♦ 1TE 

IF  (OZOXL(I.J) .LT.O.)  GO  TO  430 

1 2L ( J ) =  I 

GO  TO  450 

CONTINUE 

CONTINUE 

WRITE  ( b  *  4b9 ) 

00  460  J- 1 * JM 

WRITE  (6. 461)  J.IZU(J) .IZL(J) 

DO  500  1  =  1  * IM 
DO  500  J=1«JM 
DELU ( I  ♦  J) *0. 

DELL ( 1 1 J ) *0 • 

DO  550  1=1. IM 
DO  550  J=1.JM 
KX=XR0DU(1» J) 

IF  (KX.EU.O)  uO  TO  525 
OFLUU  tJ>  BZE  TA  (mX)  -ZETAUt  I  ♦  J) 
KX=KROUL ( 1 . J) 

IF  (KX.EO.O)  GO  TO  550 

DELL (  I . J) =ZETAL (I . J) -ZETA(KX) 

CONTINUE 

DO  600  1  =  1  . IM 

DO  600  Jsi.JM 
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IF  U)ELU(1*J)  .Gf- .0.)  RO  TO  690 

WRITE  <6*001)  I*J 

STUM 

590  IF  (DELL < i * J)  *GE  .0.  )  GO  TO  600 
WRITF  (6,002)  J,J 
STOP 

600  CONTINUE 
C 

DO  700  I  *  1  *  I M 
700  OELS ( I ) =0. 0 

00  no  1  =  IN0SE,ITAIL 
710  DELS ( I ) =GAMF 

DO  7? 0  I=iLE*iTE 
720  DEL  S ( I ) =  GAMw T 
IWX=IFS-1 
XF1=X ( I F  S» 1 ) 

ALC3=0.6*  <3.*YCAPO/ETAO-HPI*A3) 

ALC4  =  0.S*  <hPI«A3-YCAP0/ETA0) / (tTAO«ET AO) 
YFS=PF  < 1 ) 

00  722  N=2 * NF 

IF  (YFS.Ll *BF (N) )  YFS=«F(N) 

722  CONTINUE 

DO  790  I  =  INQSE  *  I FS 
J2= JROO  ( I ) 

J1=J2-1 
X1=X ( I  * J1 ) 

X2  =  X  < I  * J2) 

Y I =  Y < Jl ) 

Y2  =  Y  < J2) 

SLOPE s ( YFS-Y1 ) / ( Y2-Y1 ) 

XTtST  =  Xl ♦SLOPE* <  X2-X  1 ) 

IF  (XTEST.GF.XF1)  GO  10  *00 
XXX= ( <  X 1 -XF 1 ) / XF 1 ) **2 
YF 1=YFS*SUBT ( 1 .-XXX) 

XXX= (  ( X  2-XF 1 ) /XF  1  ) «*? 

YF2sYFS*SUB  T ( 1 .-XXX) 

PAT IOs ( Y2-Y1 ) / (YF2-YF1  ) 

YE= ( Y 1-PAT  IQ* TFl ) /( 1 ,-PaTIO) 
YCAPEsYE/t)SPAN 
E  T  A 1 =E  T  A ( U2 ) 

730  ETA1SU=ETA1*ETAI 

FUN  =  FTA1* ( ALCJ*ALC4*FT  A  ISO) -YCAPE 
FUNPR=ALCJ^3.*ALC4*ET AlSO 
DEL=FUN/FUNPR 
E  T  A2=E  T  A 1 -DEL 

IF  (DEL.LT. l.E-10)  GO  TO  750 
ETA1=ETA2 
GO  TO  730 
760  ETAE=ETA2 

GAM  =  ABS (E  TA ( J2) -ETAF) 

OELS ( I ) =GAM 
790  CONTINUE 
800  CONTINUE 
WETUPN 
ENU 


SUBROUTINE  coeff 

CALCULATES  FIXEO  GEOMETRY-  AND  GRID-RELATED  QUANTITIES  WHICH 
APPEAR  IN  ThE  TRANSFORMED  POTENTIAL  EQUATION. 


COMMON/CAPS  /  CAPG<49*’9) , CAPQT 1 ( 19*25) ,CAPGT2<49, 19) »CAPl < 19) » 

*  CAPHU9  9)  ,  CAPHTl  (19*25)  *CAPHT2(49tl9)  .CAP  10(19)  » 

*  CAPGBd  *9)  *CAPIT  (49*19) 

COMMON/COEF  /  A 1 (49 ) *  *  49 ) , A  3 (49 ) , A5 ( 49 ) *B l < 1 9) *B2 ( 19) *B3 ( 19) , 

*  B5 (19) ,B> .19) *C1 125) *C2 (25) *C3 (25) , C5 (25) »C7 (25) » 
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COMMON/CONST/ 


COMMON/COORD/ 


D1*D2#E1*E2*F1»F2 

ALPHA.ZMACH.DPLI'‘WEPSI**E*WG»bSPAN*AI2.CSA.RXl*RX2* 
SNAKDETAtTDXl .  TDZE T A . DE T A . OX  I » OZET A . iLt ♦ IM* INOSE * 

I  TAIL* ITEt JM* JROOT * JT IP »KM , <W » JB1 * JB2 » JB3 »PLCB 1* 
PLC82 ♦ PLCH3 

CHORD < 19) .OCDY  <191 *0T0Y(19) .UXLEOY ( 19)  *ETA ( 19) *F (49) 
.FLE*FN.FTE»G(19>  iGAMLE.GAMN»GAMF.GAMTE*GAM#T»H(?S>  • 


CH0wD<19) ,OCDY  <191  *UT0Y<19)  .UXLEOY  < 19) » ETA ( 19)  *F (49) 
*FLE»FN,FTE»G<19>  iGAMLE*GAMN.GAMF*GAMTE»GAM#T.H<2&>  . 
TAU(19) «X(A9*19) *XCAP(49) ,XLE(19) ♦Y(19)iZETA(251  *GF. 
ZCAP ( 25 ) *ZLE (19) »UZlEOY ( 19 ) .ALPHATU9) .0ADY<19) »GmT 


PI=4.*ATAN ( 1 . ) 

ILXsILE-1 
I  TP= I TE* 1 
nu  2nn  j  =  kjm 

CC=ChOMQ ( J) 

tt=tau ( J) 
rc=TT*CC 

<)XLY  =  OXLEUY  <  J) 

OZL  Y  =OZLEUY ( J  > 

OAT Y=OAOY <  J) *P 1/180. 

CIS=-OCOY (J)/CC 
CTS=-OTOY ( J) / fT 
CAPI < J>  =C1S/CC 
CAPIR  ( J)  =  (C1S*CTS)  /TC. 

00  100  K=1,KM 
ZZ=ZCAP (K) 

CAPGT1 ( J*K) =  ZZ* (CIS*CTS) 

I  CAPHTl ( J,K) =2.*ZZ* (CIS*CIS4CTS*CTb4CIb*CTb) 

DO  120  1 = 1 » I M 

CAPG<  I  .  J)  =  (XCAPU  >  *0.5)  *CIS-OXLY/CC 
i  CAPH(I,J)=2.*L1S*CAPG(I,J) 

ANGLE =ALPH AT ( J ) *P 1/180. 

TNAaTAW (ANGLE) 

SCA=1./C0S(ANGLE) 

00  140  1  =  1  * ILX 
CAPGR ( I » J ) =0 • 

CAPGT2<I*J)=-0ZLY/TC 

CAPHT2<I*vJ)=2.*CAPGT2(I  ♦  J)  *  (CIS*C  IS) 

S^UoMalLEilTE 

XX*XCAP(I)  ♦0.«5 
CAPGB(1*J)*TNA/TC 

CAPGT2< I . J)=-(OZLY-CC*XX*SCA*SCA*DATY*DXLY*TNA)/TC 
CAPHT2( I  * J) =-2.*SCA*SCA*0ATY*(0XLY-CC*XX*fNA*DATY) /TC 
1  ♦2.*CAPGT2(I»J)*(CIG*CTS> 

i  CAP  IT  ( I.  J)  =SCA*SCA«DATY/TC*CAPGt5(  I  » J)  «  (CIS*CTS) 

00  180  1  =  1  TP* Im 
CAPGR ( I » J) =0. 

CAPGT2 ( 1  * J) =- (OZLY-CC*SCA*SCA*UATY*OXLY*TNA)/TC 
CAPHT2 ( I » J) =2»*SCA*SCA*DATY* (DCUY ( J) *CC*TNA*DaTY)/TC 
1  ♦2.*CAPGT2 ( I , J) • (CIS*CTS) 

i  CAP  I T ( I « J ) =0  • 
i  CONTINUE 

OX*2.*DXI*OXI 
0E=2.*0ETA*0ETA 
OZ=2.*OZ£IA*OZETA 
A l (  I ) =0  • 

AKIM)  =0. 

A2<1)=0, 

A2 (  IM) =0. 

A3( i) =0. 

A3 ( IM>  =0. 

IMA=IM-1 

00  2 l Q  I=2»IM* 

A1 (I)  =  (F ( I  ♦ 1 ) ♦ F (I) ) /OX 
A3 ( I ) * (F ( i ) ♦F ( I  —  1 ) ) /OX 
A2< I ) s-Al ( I ) -A3  (  I  ) 
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AS ( 1 ) =0. 

AS  ( 2 )  =0. 

00  £20  1  =  3,  IM 

?2n  A5<I)  =  (F{i-l)*F(W> )/OX 
C 

Si  (  1 )  =  <G ( 1 ) *G (2) ) /OF 
S3  (  1  )  =  (G(  1  )  *0(2)  )  /he 
S2(l ) =-Hl ( 1 ) -S3 ( 1 ) 

31  ( JM) =0. 

32  < JM) =0. 

B3(JM>=0. 

JMX=JM-1 

1)0  310  J=2»JMa 
31 ( J )  =  ( G ( J*  1 ) ♦G(J) ) /nt 
S3  (  J)  =  <G(  J)  *C,  (  J-l  )  )  /HE 
310  32<J)=-ril (J)-83(J) 
H5(1)=(G(2)>G<3) )/OE 
35 (2)  =  <G<  1 ) *G(2) ) /OE 
00  320  J  =  3  ♦ JM 

320  S5(J)=(G(J-l>*G(J-2>)/DE 
JMX2= JM-2 
00  330  J=1»JMX2 
330  37(J)=(G(J*1)*G(J*2) )/OE 
37 ( JM-1 ) =U. 

S7 ( JM) =0. 

c 

Cl  (l)=u. 

Cl (KM) =0. 

C2( 1 ) =0. 

C2<KM)=0. 

C3 ( 1 ) =0  . 

C3  <KM) =0. 
kmx=km-i 

00  410  K=2.KMA 
Cl (K) = (H (K* 1 ) *H (K ) ) /OZ 
C3(K)=(H(M *H(K- 1 ) >/0Z 
410  C2  (  K )  =-C 1  (K)-C3(K) 
CS(1)=0. 

C5  <2) =0. 

00  420  k  =  3  »  KM 

420  C5 (K) = (H (K-l ) *H (K-2) ) /OZ 
KMX2=KM-2 
00  430  K=i»KMX2 
430  C7 (K) = (H<K*1 ) *H (K*2) ) /DZ 
C7 (KM-1  ) =U  . 

C7 (KM) =0  . 

C 

01  =  1  ./  <4.*0X1<>DETA) 
02=-01 

El=l ./ (4.*0ETA*0ZETA) 
E2=-F1 

F  1  =  1 .  / ( 4 . *0X I  *0ZET  A) 

F2  =  -F  1 
return 
END 


C 

c 

c 

c 

c 


SUBROUTINE  INlT(MSTARf) 
INITIALIZES  solution  ARRAYS. 


» 


COMMON/CONST/ 


COMMON/SOLVO/ 


ALPHA.  7MACH»r>PLI" 
SNA. TOE  TA.TOXI. Tl 


iM.EPSI .WE.WG.BSPAN. AI2.CSA.HX1 .RX2» 
.  _  .  luZE TA.DETA.OXI , OZE T A » I LE » I M ♦ I  NOSE » 

I  TAIL. TTE . JM, JROOT. JT IP.KM.KW. JS1 . JB2. JB3.PLCB1 . 
PLCH2.PLCB3 

ULl  (49.  IB)  ,0L2(4  9, 16)  .DSl  (4V)  ,0S2(4R)  .01)1  <49»16)  f 
0U2 (49.16)  .QPMAX.OPSAXT .GAMMA (16) , I MAX  I , I  TEST , 

JMAXI , JM A XT, KM AX  I » K M A  X T ♦ NSUP ♦ P < 4 9 , 1 6 . 25 ) , PLE1  ( 16) t 
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00,0 


*  PNEw(PS)  ,PN0SE,PSYM(49,P5>  ,PT1  $16.25)  *PT2<16,25>  » 

*  PT3U  6,25)  »PTE1  (16)  *PTE2  (  lb)  »PwRL  (49.  i  6)  ,PwBU(49»16> 
C0MM0N/5UHF  /  DELL (49*16) * DLLS (49) *UELU(4** 16) * 02DXL ( 49 « 1 6 ) * 

*  f)Z0xU(49»  16)  ,DZ0YL(49*16)  * OZUYU < 49 , 1 6 )  ,rlWBL<4  9,16)  * 

*  HwBtj  (49*  16)  *1ZL(19)  *I2U(16)  ,UHOO<49)  »KB00L  (49*16)  » 

*  KR0QI1 ( 49  » 1 6 ) *ZL (49,16) *ZU<49»16) 


100 

106 

110 


120 

uo 

140 


145 

150 

300 


700 

750 

800 

900 


IF  (MSURT.NE.O)  60  TO  105 
00  100  1  =  1  «1M 
00  100  0=1 *JM 
DO  100  6=1, KH 
P<  1*  J.K)  ai), 

DO  110  J=1*JM 
DO  110  6=1»KM 
Pfl ( J*K) sU, 

PT2 ( J*K) =0. 

PT3(J*K) =0. 

00  IPO  1=1. IM 
DO  120  J=1»JM 
0L1 (I.J)sO. 

0L2( I ♦ J) =0. 

PwBL ( 1 » J) =0. 

DU  1  ( I » J) =0. 

OOP ( I « J) =0, 

PwBU ( I ♦ J) =0. 

DO  130  1  =  1,  IM 
DS1 (I)=0. 

052 ( I ) =0 . 

DO  140  J=1,JM 
PLE1 (j;=0. 

PTE  1 < J) =0. 

PTE2 ( J) =0. 

PNOSF=0. 

00  145  1  =  1,  IM 
00  145  K=1.KM 
PSYM  < I ,K) =0. 

IF  <MSrART.NE.O)  GO  TO  300 
00  150  j=l,JM 
GAMMA ( J) =0. 

P I =4 . * AT  AM ( 1  •  ) 
AMGLF=ALPMA*PI/1«0. 
CSA=COS< ANGLE) 

SMA=S IN (ANGLE) 

PX1=1 ,/wt 

PX2=1 .-HXl 

T()XI  =  2.*0A1 

TOE  T  A=2  »  4'  'ETA 

TDZETA  =  2,*QZEU 

A I  2=1 ./ (ZMACH*ZMACH) 

IF  (MSTARl .EO.O)  GO  TO  900 
00  HOO  1  =  1 NOSE  *  I M 
JEX=JHOO( I )*1 
DO  760  J= 1 » JEx 
KL  =  KRODL ( I »  J ) ♦ 1 
KU=KROOU( 1 ,0) -1 
00  700  k=*L»KU 
P( I , J.K) sO, 

CONTINUE 
CONT  FNDE 
RETURN 
ENO 


SUBROUTINE  SOLVE (ITER) 

£  EXECUTES  A  COMPLETE  RELAXATION  SWEEP  OVER  THE  COMPUTATION  DOMAIN, 
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nnonri 


i 

C 


DIMENSION  Sl**l  <49. 19) ,SL2( 19) *SS2<19) *SU2<19) 
DIMENSION  SUK#1  ( 49 . 19) *SL3 < 19) »SSi ( 19)  *SU3(19) 


COMMON/CAMS  /  CAPG < 49. 1 9 ) * C APGT 1 ( 19. 25) .CAPGT2 < 49* 1 9) .CAP  I < 19) ♦ 

C APH (49*19) .CAPHT1 U9.25) .CAPHT2 (49* 1 9) ,CAPIB(19>  * 
CAP6B (49*19) ,CAP1? (49.19) 

ALPHA,ZMACH*DPLIM*EPSI*WE*W5»8SPAN,A12»CSA*PX1*RX2* 

SNA  .tdeta.tdxi.tozeta.deta.oxi.ozeta'Ile.im. inose  * 


COMMON/CONST/ 


I TA IL  » I 7£ ♦ JM* JROOT  .JTIP.KM.KW.JBl » JB2*  J83.PLCB1 « 
PLCB2*PtCB3 

COMMON/COORO/  CHORD (19) ,DCDY<19) *DTDY(19) »DXLEDY(19) »ETA(19) *E(49) 


COMMON/SOLVO/ 


»FtE*(rN.PrE*G(I9)  *GAMLE*GAMN»GAMr*GAMTE*GAMwT.H(2^)  . 
TAU(19)»X(49.19) *XCAPT49) *XLE(19) ♦ Y { 19) ,2ET A ( 25) .GF ♦ 
ZCAP  ( 25)  »  ZLE  (19)  .0ZlEDY<19)  ♦AIPHATH9)  ,0A0Y(19)  .G#T 
DL1 (49.19) ,DL2(49» 19) .QSl (49) *DS2(49>  *DU1 (49.19) ♦ 
DU2(49.19) .DPMAX.DPMaXT*GAMMA(19T*IMAXT ,lTERt* 

UMAX I »  JMAX  T .KMAXl ♦ KMAXT . NSUP  *P (49*19*  2b ) »PL£1 (19) * 
PNEtf (25) .PNOSE .PSYM ( 49. 25) .PT1 (19.25) .PT2(19,25) . 
PT3( 19.25) *PTE1 (19) .PTE2( 19) »Pw8L(49*l9>  *P*BU (49*19) 
COMMON/SURF  /  DEI  L(A9,19) .DELSIA9) .0ELU(49.19)*DZDXL(49*19)* 

DZOXU(4  9, 19) .DZDYL (49.19) *DZDYU<49* 19) .H W8L(49,19)  ♦ 

!S8)Ut2«!nr  !2fc  »*»  • 


UPDATE  SURFACE  CONDITION 

UPMAX=0.0 

NSUPsO 

ilx=ile-i 

I MX  =  I M-  1 

INXalNOSE-1 

ITP=ITE*1 

JMX=JM-1 

KHX  =  KM-  1 

DH=-OX  I 

OG=-GAMN 

Cl=2.«OH«UH/ (UH*2.*nG) 

C2=4. *DG/ <DH*2, *DG) 

C3= (OH-2, *DG) / <DH*2,*UG) 

PW  =  P< lNUSt-1 *  1  *KW> 

PWW  ;P  ( I  NOSE -2*  1  *K«l) 

PNOSE  =Cl*CHORL'  ( 1 )  *CS&/FN*C?«PW*C3*PWW 
OGs-GAMLE 

Cl*2.*UH*DH/ (uh*?,*OG) 

C2  =  <*» *OG/ <0h>2. *DG) 

C3*(OH-2,*OG)/ (DH*?.»OG) 

CF=CS A/FLt 

00  20  JsJROOT.JTIP 

Pw*P ( ILE-1 » J*NW) 

PWw=P  (  ICt''2*  J*KW) 

20  PLE1 ( J) =C1*CH0R0(U) *CF*C?»P**C3*PtfW 

PLE1 ( JROOT-1 ) =2.*PLE1 (JROOT)-PLtl (UkOOT*l> 

00  30  MSUkF s  1  *  2 
DO  30  ISFIsINOSF.IM 
JE  =  JHOO ( I  SET ) *  1 
00  30  JSE  T* 1 «  JE 

30  CALL  SURFbC ( 1SET . JSFT .MSUhF) 

00  4  0  I SE  f  s I NOSF *  I M 
JSETsJrfOOlISET) 

40  CALL  SUHE'SC  ( 1  SET  « JSFT  *  0  ) 

CALL  S YmRC (0*0*0) 

OHsOX  l 
))G*GAM  T  t 

C1*2.#0H*UH/(UG*  (()H*OG)  ) 
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ono  n,  nor. 


C2=-2.M0H-DG> /DG 
C3=  (Oh-uG)  /  (OH ♦Of*) 
oo  so  o=jKoor»jrrp 

BTE=(2.-0G/0H)  *PWbL  U  TF-*  J)  -U  ,-OG/DH)  *PWHL  ( JTK-1 »J) 
pe=p ( i  re*  1 1 jtKw) 

PEE  =P  ( ITE*2. J»Kw> 

PTE  1 ( J> =Cl*PTt*C?*PF*C3*PfF 
so  PTfc2(j)=3.*(PTei<j)-P£i*pee 

PTE  1 ( JR0QT-l)=2.*PTEl (JROUT)-PTEl (JR00T*1> 

COMPUTE  POTENTIAL:  REGION  UPSTREAM  OF  FUSELAGE  NOSE 

P ( INOSE ♦ 1 *KW) =PNOSE 
IF  UI2.GT.1.)  GO  TO  103 
00  102  J=1»JM 
HO  102  K=1«KM 
PT3(J»K>=P(1.J.K> 

102  PT2(J.K)sP(2*J«K) 

ISTART=3 

GO  TO  10S 

103  00  104  J=1.JM 
00  10 4  K=1*KM 

104  P  T2 ( J  «  K ) =P (  1  ♦ J,K) 

1ST  APT  =  2 

105  00  ISO  IStT=lSTART« INK 
00  120  JSET=1«JmX 

CALL  TRICOE (2«KMX,1SET«JSFT) 

CALL  INVERT <2«KMX) 

CALL  MAa1<2.KMX.  ISFT*  JSET) 

00  US  Ks2 f KMA 

PT1  ( JSE  T . K ) =P (  ISFT  *  USE  T  »  K  > 

US  P  ( ISE  T.  JStT.K)  =PNFW  <K) 

PT1  <JSETU)=P<ISET«JSETU) 

PT1 ( JSET.KM) ( ISET. JSET.KM) 

120  CONTINUE  , 

00  125  K  a  1  ♦ KM 
1  25  PT 1  ( JM ♦ K ) =P ( ISET ♦ JM  »K ) 

00  130  J=ltjM 
00  130  K=1.KM 
PT3<J.K)=PT2(U*K) 

130  PT2(J»K)=PTl  (J«K> 

ISO  CONTINUE 

COMPUTE  POTENTIAL:  REGION  DOWNSTREAM  OF  FUSELAGE  NOSE 


00  200  J= 1  *  JT  IP 
200  P ( ILE.J.KW) sPLElt J) 

00  210  Jsl.JM 
SU2  < J>  =PT2 ( J.Kw) 

SU3 (J) =PT3  < J.KW) 
SL2<J>=PT2<J.niv> 

SL3U)  *PTJ(J.Kw) 

SS2  ( J )  =PT2 IJ.KW) 

210  SS3(J)=PTJ(J.Kw) 

PSAVE  =  P  ( ILE.  JH00T-1  iN) 


00  bOO  1  St  T  =  INOSF  *  I  MX 

isetx=isf.t-i 

I SE  TP= I SF  T ♦ 1 

JE2= J800 ( I SE  TA  > -  2 

JE 1 = J800 ( 1  St  T ) -  1 

UFO=JPOO(ISETP)-1 

IF  (ISEI.tQ.lLX)  JfOaJEl 

IF  (ISET.EQ.I TP)  j6  2  =  UE  1 


LINE  SEGMtNTS  BELOW  WING/HOUT 

00  220  JsltJM 
PT2(J.KW)=SL2(J> 

220  PT3 ( J.Kw) =SL3 ( J) 

IF  (JE2.E-..0)  GO  TO  230 
DO  225  J=i * JE2 
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c 

c 

c 


SUKW]  <  ISE  TX,  J)  =  P  <  ISF  TX,  J,KW*  1  ) 
K8=KR0OL  (  ISF.Ta,  j) 

P(  ISFTX,j,KB*i  )  =OU  (ISETX,J) 

22S  P  ( ISETX,  J,KB*2)  =')L2  ( ISETX,  J) 

2 JO  CONTINUE 

DO  235  J= i * JE 1 

SUKW1 (ISEf,J)=P(ISET»J»KW,l) 

KB  =  KRODL ( 1 SE  T *  J ) 

PIISET, J,Kb*l»=DLl (ISET,J> 

2J5  PUSETf J,KB*2)=Dl2(ISET,J> 

00  240  J= 1  * JEO 

SUK*  1  <  ISE  IP*  J)  =P  (ISF  TP*  J,KW*  1  ) 
<B=KHODL ( ISETP, J) 

P  ( ISETP.  J,KB*  i )  =I)L1  ( ISE  TP,  J) 

240  P(lSETP.J*KB*2)=nL2(ISETP.J> 

DO  260  JStT  =  1 » JE  1 
Klj  =  KRODl  (  ISE  T.JSFT) 

CALL  TRICUE  (2,KH, ISET, JSFT) 

CALL  INVENT ( 2  .KB) 

CALL  MAXI (2. KB, ISET, JSET) 

DO  245  K=2,KR 

PT1 ( JSET  ,K) =P (I SET, JSET,K) 

246  P(ISFT»  JStT»K)=PNF.«(K> 

PT1  (JSET,  i)=P(  ISF.T  ,  JSET.  1) 

K8 1 =KB*  1 
00  247  K=KB1,kw 

247  PT1 (JSET, K)=P ( ISFT, JSET, K) 

CALL  SURFdC ( ISET, JSET,  1 ) 

P ( ISET, JStT » Kri ♦ 1 ) =0L 1 ( I  SET , JSET) 
P ( I  SET, JStT ,K»*2) =DL2 ( I  SET , JSET) 
250  CONTINUE 

DO  265  J= 1  *  JE 1 
SL3 ( J) =SL2 ( J) 

255  SL2  <  J)  =PT1  ( J ,m») 

LINE  SEGMENTS  ABOVE  wlNG/HODY 


DO  260  J=1,JM 
PT2 ( J.Kw) =SU2 ( J) 

PT3 ( J , Km ) =  SU  J (J) 

IF  (JE2.E'-.0)  GO  TO  210 
DO  265  J= 1 , JE2 

P(ISETX,j,NW*l)  sSUKWl  ( ISF.  TX,  J) 
SLKwl  <  ISE IX, J) =P ( ISETX, J,KW-1 ) 
KB=KBODU ( ISETX, J) 

P  (  ISETX, J, KB- 1 ) =0U1  ( ISETX, J) 

P  ( ISETX, J.KB-2)  =0U2 ( ISETX, J) 

CONTINUE 

DO  275  J= i , JE i 

P(  ISFT »  J.KW, 1 ) =SUKw 1 < ISET , J) 
SL**1 ( ISET, J) sP ( ISET, J,KW-1 ) 
KB=KROJU( I SE  T , J ) 
P(ISFT,J,KB-i)=DUl ( ISET , J) 
P(ISFT,J,KB-2) =DU2 ( I SE  T , J ) 

00  280  J= I , JEO 

P(  ISETP, J.Kw* 1)=SUKW]  (  I SETP , J ) 
SLKwl <ISE TP, J)=P(ISETP,J,KW-I ) 
KB=KBDDU ( ISETP, J) 

P  (  ISF  TP,  J.KB-l )  =T'U  I  ( ISETP,  J) 

P(  ISF TP, J,K6-2) =0U2< ISETP, J) 

00  290  JSE T  =  ’  , JE  1 
K8=KP0QU ( ISET, JSET) 

CALL  TRICuE (Kp, KMX, ISET, JSET ) 
CALL  INVERT (KB , KMX ) 

CALL  MAXI (KR, KMX, ISET, JSET) 

DO  285  K=KB,KMx 

PT1  ( JSE  T , K ) =P ( ISFT, JSE  T , K ) 

P  (  ISFT, JStT ,K) =PNE 4 (K ) 

PT1  (JSE  I  ,KM) sP ( ISFT , JSET, KM) 

KH1=KR-I 

DO  28  7  KrKW.Kdl 
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?ri  7  PT1  <J5FT,K)=P  (  ISET, JSET. K) 

CALL  SUPFdCUSF  T.JSkT.2) 

P ( I SF  T , JSt  T  « Kri-  1 ) =  0U1  ( I  SET  , JSET ) 
P  ( 1ST  T*  JStT  ,Kd-2>  =  DU?  ( ISET.JSED 
29ft  CONTINUE 

i)0  295  J= l , JE  1 
SU3< J)  =  SU2  t  J) 

29S  SU2  ( J )  =PT  1  <J,Kw) 

IF  (JE2.EvJ.0>  GO  TO  310 
DO  30b  J=  1 *  JE2 

30S  P<  ISFTX,  J.KW-l )  =SLKt»l  (ISETX.J) 
310  CONTINUE 

DO  315  J= 1 »  JE  i 

3 1 S  P( ISET. J,Kw-l ) =SLKW1 ( ISET, J) 

DO  320  J=i«JEU 

320  PdSFTP.JtKw-naSLKwl  (ISFTP*J) 


LTNES  0 U T d 0 a P u  OF  WING/RODY  SlOt  EDGE 


325 
332S 

3330 

326 


3327 

327 

330 


33? 


335 


35ft 

370 
3  7S 

400 

41ft 


P< ILE. JWOOT-1 *KW) =P5AVE 
GO  TO  332S 


IF  (ISET.tO.lLX) 

IF  ( ISET  ,tO.  I  TP) 

DO  325  J=1»JM 
PT2(J,Kw) =SS2 ( J ) 

PT3< J.Kw) *SS3<J) 

GO  TO  326 

DU  3330  J=JROOT.JTIP 
PT2(J.K*)=PTE1  (J) 

PT3 ( J.Kw) =PTE2 ( J) 

P( ITE, J»KW)=PTE1  <J) 

P ( I TF-1 » J*KW) SPTF2 ( J) 

P ( I TE» JHOUT-1 ,KW) =PTF 1 ( JPOOT-1 ) 
IF  (JE2.EO.O)  GO  TO  327 
IF  ( ISET.tU. I  TP)  GO  TO  3327 
P ( ISFTX, Jt2.*w) =DS1 ( ISETX) 

IF  (JE2-1.LT. 1)  GO  TO  327 
P ( ISFTX, Jt2-1 ,K*) =DS?( ISETX) 


CONTINUE 

PT1 (JFl.K4)=nSl (ISET) 

IF  (JE1.LT.1)  GO  TO  330 
PTl(JFl-l»Kw)=US2(ISFT) 
CONTINUE 

IF  ( ISET.EU. ILX)  GO  TO  332 
P< ISFTP. JtO.KvO  =0S1 (IbETP) 

IF  (JEO-l.LT.l)  GO  TO  332 
P< ISFTP,JtO-l »Kw) =QS2< ISETP) 
CONTINUE 
JSET=JBUf)(ISEf ) 

CALL  SUHFSC ( ISETt jSFTtO) 
P(ISFT.JE1.KW)=US1  (ISET) 

IF  (JE1-1.LT.1)  GO  TO  335 
P( ISET. JEA-1 ,Kw) sDS? (ISET) 
JE= JF 1*1 

00  370  JSt T  =  JE  « JMX 

CALL  TRICUF (2,KMX» ISET* JSET) 

CALL  INVERT (2»KMX) 

CALL  MAXI <2, KMX, ISET, JSET) 

no  350  K  =  2 ♦ KMX 

PT1 < JSET ,K) =P ( ISET. JSET, K) 

P ( ISET, JSET, K) =PNEW (K) 

PT1 <JSET,1)=P< ISET.JSET.l) 

PT1 (JSET ,KM) =P ( ISET ,JSET ,KM) 

CONTINUE 

00  375  K  =  1 , KM 

PT1 < JM,K)=P ( ISET, JM,K) 

00  400  J=1,JM 
00  400  Kri, KM 
PT3(J,K)=PT2(J,K) 

PT2 ( J,K) =PT 1  (  J,K  ) 

00  410  J=1»JM 
SS3 ( J ) =PT3( J.KW) 

SS2  <J)=PT2(J,Kw) 
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00  420  J  =  j£  i 
$L? ( J )  =SS2 ( J) 

420  su? < j) =552 < j) 

500  CONTJNUt 

wRITf  (6.511)  I  TFR.DPMAX. 1MAX I  * JMAXI .KMAX1  ,NSUP 
511  FORMA r  (in  ,110, E15. 4,  3 15. 110) 

UPDATE  CIRCULATION  AN'U  FAW-FIELU  CONDITION 

DO  oOO  J=JHOOT . JT1P 

f,NFw  =  (2.-OG/DH)  *  (PWHU  (  I  TF  .  J) -PWttL  (  ITE  .  J)  )  “(  1  ,-DG/DH) 
1  (PWHU  l  I  Tt- 1 «U ) -PwHL ( I TE- 1 . J)  ) 

600  GAMMA ( J) =wG*GNEw* ( 1 ,-«6) *GAMMA (U) 

GAMMA ( JPOUT-1 ) =2.*GAMMA ( JPOOT) -GAMMA ( JROOT  *1 ) 

CALL  FAwRC 

wPITf-  (6,oll)  GAMMA  (  JMOOT)  *  I  TEH  1  .DPMAXT  *  JMAXT.KMAXT 
611  FOWMAT  (M*.  55  X.E15. 4.110,  E  13.4.215) 

WE  TOWN 
END 


SUBROUTINE  FAW0C 

UPDATES  FAR-FIELD  BOUNDARY  CONDITION  IN  THE  TREFFTZ  PLANE. 


dimension  save ( i 9 ) 
COMMON/CONST/  ALPP 


ALPHA.ZMACH,OPLlM»EPSI .WE.WG.BSPAN.AI2.CSA.RX1.RX2. 
SNA.TOETA.TOXI .TUZETA.DETA.JXI.UZETA. ILE.IM.INOSE. 
ITAIL.ITE. JM. JROOT. JTIP.KM.Kw.JB1.JB2. JB3.PLCB1. 
PLC82.PLCB3 


COMMON/SOLVO/ 


KM.Kw, JB1. JB2. JB3.PL( 


DL1 (49.19) ,DL2 (49.19) ,DS1 (A9) .DS2(49) .DUl (49.19) ♦ 

0U2 (49. 19) .DPMAX .DPMAXT .GAMMA ( 19)  . I  MAX I *  ITERT  t 
JMAXI , JMAXT.KMAXI .KMAXT.NSUPtP(49. 19.25) »PLE1 ( 19) ♦ 
PNEW ( 25) .PNOSE.PSYM (49.25) .PT1 (19.25) .PT2( 19,25) » 

PT3 (  19,25)  , PTE  1 (19) tPTE2(19) .PwBL(49,19) ,PW3U (49,19) 
DELL (49»19)  .DELS (49) ,DELU( 4^.19) .DZDXL  U9, 19) » 

DZDXU (49, 19) .DZDYL (49, 19) .DZDYU (49. 19) .HWBL (49, 19) , 
HWBU ( 49 .19) » I ZL ( 19) .IZU<19> »JB0D(49) .KBODL < 49, 19) * 
K80DU  (  49 »  19)  .ZH'49.19)  .ZU(49,19) 


I TEWT=0 

JMX=JM- l 

K  MX  =K  M- i 
JE  =  JRUD ( I M ) - 1 
100  OPMAXT  =  0.0 

00  110  J=l.Jt 

SAVE ( J) =P ( IM, J.KW.l ) 

KR=KROOL ( I M , J ) 

P ( I M , J ♦ KR ♦ L ) =OL 1 ( IM. J) 

110  P ( IM, J,KB*2> SDL2( IM. J) 

DO  150  JsI.JF 
KH=KB0DL ( IM, J) 

JSE  T  = J 

CALL  TWIT  (2.Kri,0.jSET) 

CALL  INVERT*?, KP) 

00  125  K  =  2 , K B 

UFLP=PNEw(K)-P(IM,j,K) 

()P=AP5  (DELP) 

IF  ( f)P , L T  •  UPM u x T  )  GO  10  125 
QPMAXT  =  op 
JMA  X  T  =  J 
KMAXT=K 

PT1 ( J.K) =P ( JM, J.K) 


75 


nnnnn 


9 


US 


13(7 

150 


210 


225 


230 

250 


310 


P ( IM, J»K) =PN£w  (k  ) 

PT 1 < J, 1 >  =P l]M, J, 1 > 

KB1=KB*1 

(70  i  30  K=*81»>\w 

P71 (JiK) =P (  2 M »  J » K ) 

CONTINUE 

00  210 

P  (  IM»J*KW*1 )  =  5  A  V  E  ( U ) 

SAVE ( J) «P l lM»U»Kw-l ) 

KRzKBODU ( 1M,J> 

P ( IM, J,KR-1 ) =UU1 ( IM, J) 

P  (  I  M  ,  J  *  Kfc)-2  )  =0U?M  I  M  ,  U  ) 

00  250  J=1*JE 

KH=KR00U ( IM, J) 

JSETsJ 

CALL  T«I  T  <KB,*MX,0,  JSt'T) 

CALL  INVtHT  (KB,KMX) 

00  225  K=K8,KMx 
0£LF=PN£w (K|-P(JM ,j.k; 
i)P“  APS  (  OF  LP  ) 

IF  (OP.LT .OPMAXT I  GO  TO  225 
0PMAXT*0P 
UMAX T  =  J 
KM A  X  7  =  K 

PT1  < J«K) aP( [M, U,<> 
P(1M,J,K)=PnE*(K) 

PT 1 < J,KM) =P ( IM, J,KM) 

KB1=KB-1 
DO  230  KaKWtKB 1 
PT1  ( J,K) =  P ( IM, J.K) 

CONTINUE 

DO  310  J=  1  » JE 
P<IM,J,K*-1)=SAVE(U) 
P(IM,UE,KW)=DS1 (IM) 

IF  (JE.NE.l)  P(IM,JE-1,KW)sOS?(IM) 
PT1 (jE.KW)sDSI (IM) 

IF  (JE.NE.l)  PT1(JE-1*KW)=DS2(IM) 
DO  350  U=UH00T ,JMX 
JSt  t  = j 

CALL  TBIT <2, KMX, 0, US* T) 

CALL  INVENT (2, KMX) 

00  325  K=2,KM* 

DELP=PNEW (K) -P ( IM, J,K) 

DP=APS (DELP) 

IF  (QP.LT.OPMAXT)  GO  10  325 
DPMAXT  =  DP 
JM AX  T  a J 
KMAXT  =K 

PT1 ( J,K) =P ( IM, J,K) 

P(IM,J*K)=PNE«*  { K  ) 

PT 1 (U*1)=P<IM,U,1) 

PT1 < J,KM) =P ( IM, J,KM) 

CONTINUE 

ITEHT=1TEHT*1 

IF  (  IT£«T,ut  ,^00>  RETURN 

IF  (nPMAXT.oT.UPLIM)  GO  TO  100 

pe turn 

ENO 


SUBROUTINE  SURFBC ( I , JtM) 

CALCULATES  INTERIOR  IMAGE-POINT  POTENTIAL  VALUES  REQUIRED  TO 

satisfy  the  flow-tangency  bounoa«y  condition  on  lore*  <m.eq.i> 

AND  UPPER  <M.fcQ,2)  SURFACES  OF  THE  WING/BODY, 


76 


f 


c 

c 


8 
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COMMON/C ArS  /  CAPG (49, 1 9)  ,CAP6T1  ( 19* 2b)  , C APGT ?  ( 49 , 1  9 )  ♦CAPI  (19)  * 

►  c5pGBT29!?i?5cSpitjil:f|| -CAPHiir^9,lc,)  ,CAPishl»^ 

:C0M"0N/C0NST/ 

‘  PLC82!plCB3M*  JR°0T  JB1  » JB2, JB3»PLCB1 , 

CH0HD]l9)  ,DCDY  <  1 9)  »0T0Y<19)  *OXL£DY  (19)  #£  T A  ( 1 9)  ,F  (49) 
jFLE  jFN*FTE»G(  19)  ♦(jAMLEtGAMN*GAMF»GAMTE»GAMwT*H  (?5)  » 
*,x,9a£,(4?)  «*»-E(19)  *  Y  (19)  ,Z£TA(25) »GF, 
2CAP(25>  * ZLE  <191 »UZl£UY<19) *ALPHAT(19) »DADY ( IV) *GwT 
OL 1(49*19) «  DL2 ( 49  ♦  19 )  » OS 1(49) ,QS2 (49) ,QU1 (4  9*19)  » 

0U2 ( 49  *  19) «OPMAXjOPMaXT* GAMMA ( 1 9)  *  I  MAXI » I TERf  » 

JMAXI* JMAXT»KMAXl »KMAXT»NSUP*P (49*19*25) *PLE1 ( 19) » 

HWBU ( 49 « 19) *I2L(19)*IZU(19) *JBOO(49) *K800L (49* 1 9) * 
KBOOU (49,19) *ZL (49, 19) *ZU (49* 19) 


common/coowd/ 


COMMOM/SOLVO/ 


COMMOM/SUWF  / 


[)K  =  OE7A 

rc=T&u  < j) *chopo ( j) 

IF  (M.EG.2)  GO  TO  100 

KP  =  -l 

UL=“OZE T  A 

OEL  =  -DELL < I  * J) 

I  Z= I ZL ( J) 

SLOPX  =  OZOAL  (  I  *  J) 

SLOPY=UZDYL  (  I  *  J) 

*S  =  KRODL ( 1  * J) 

HWB=H*BL  (  i  *  J) 

CAPGTH= (OL-OEL) o (CAPGT 1  ( J » Kb > -C APGT 1  ( J,KS*1 ) ) /DL  *04861 1 ( J,KS* 1 ) 
CAPGTe=CAPGT8*CAPGT2 ( I , J) 

GO  TO  120 
100  KP=*1 

ol=dzeta 

()EL  =  0ELU  ( 1  ,  J) 

1Z*IZU( J) 

SLUPX  =  OZO*U  < I » J) 

SLOPY  =  OZOYU ( I  * J) 

KS  =  KBOOO( I  * J) 

Hv*6  =  H**B(J(  I  »J) 

CAPGTB= (DL-DEL) * ( C APGT 1 ( J ,KS ) -CAPGTl (J,KS-1 ) ) /DL*CAPGT 1 ( J*KS-1> 
CAP6T«  =  CAHGTB*CAPGT2 ( 1  * J) 

120  IF  (I.EG.1M)  GO  TO  210 

IF  U.GE.JROOn  GO  TO  UO 
QH  =  - OX  I 
IP=-1  . 

bO  TO  200 

1 4 0  IF  U.GE.1Z)  GO  TO  160 
OH  =  -()X  I 

IP  =  -1 

GO  TO  200 

160  DH=*OXl 
I  P=  ♦  1 

20  0  CAPK= (SLOPX/CHUWfX  J) ♦CAPG1 I *0) *5LOPY ) ( 1 ) 

C  A PL  = SLOP  Y  *6 ( 0 ) /RSPAN 

c  APR  =  HiVb*  (CaPGR  (  I  ,  J>  *SLOPX*CAPG(H*SLOPY-1  ./Tc) 
CAPP=CSA*bLOPX-SNA 
GO  TO  21b 
21 o  C6PK=0.0 

CAPL=SLOPY*G(-J)  /HSPAN 
CAPMsHWB* (C&PGTH*SLOPY-l ./TO) 

CAPP=-SOA 
GO  T()  220 
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21S  A 1  =- 1 . 5/i)n 
A2=2.0/DH 
A3=-0.b/0H 

uO  ro  22b 

220  A  1 =  0 . 0 
A2  =  0.0 
A3=0.0 

22b  ai=-l .b/DK 
d2=2.0/UK 
d3=“0.S/0* 

Cl=-  (nL*2.*0EL)  /  <[)tl«  (OLM)EL)  ) 
C2=<r)L*DEL)/ <l>FL*0L> 

C3  =  “0EL/(0LMUL*DfL)  ) 

01=U.S*  lDt-*0EL)  *  (?.«DL*0FL>  /  <DL#0L) 
0?  =  -nF.L*  <2.«0L*nEU  /  ( 0L*0L  ) 
03=u.b*UEL*  (UL*0FL>  /  <  L)L*DL  ) 
£1=2.*0L*UL/  (Ut.L*  (DL-MJED  ) 

E2=-2.*  (DL-DfcL  )  /OFl 
E3  =  (OL-OEU  /  <Ul*OF.L) 


Jb  =  J*  1 
JG  =  J*  2 
KS 1 =KS**P 
Kb2=KS 1 ♦  KP 
P1=P  < I »  J  *  KS ) 

P2  =  P( I  * J.0S1 ) 

P5=01*P ( I > JS»*S) *D2*P 1 1 » Jb»KSl ) *0  3*?  <1 » J5*KS2) 
P6=01«P  ( I » Jb«*S)  ♦  0 2 P  <  I  ♦  J6.KS1 >  *03*P  ( I  * J6»KS2) 
IF  IJ.E0.GH1-1.QP. J.EU.JH1 >  GO  TO  181 
IF  (J.EQ. JB2-1.0R. J.EQ. JB2)  GO  10  1182 
IF  < J.EU. JB3-1 .OP. J.EO. J83)  GO  TO  182 
GO  TO  118b 
181  JK=JP1 

PLOPLCH1 
GO  TO  183 
118?  JKzJHi 

PLt=PLCb2 
GO  TO  183 
10?  JK=JB3 

PLC=PLCB3 
183  CONTINUE 

PA=<2.*P ( I  * JK* 1 ,KS) ♦ ( 1 .-PLC) *(P<I » JK-1«KS> 

1  -?.*P ( 1 . J8.K5) ) ) / ( 1 .♦PLC) 

PB=  <2.“P  <  1 »  JKM  .Ksl  )  ♦  (1  ,-PLC)  *1P  1 1 » JK-1  »KS1) 

1  -2.*P(ItJR.KSl)))/n.*PLC) 

PC= <?.*P< 1 ♦ J«* 1 »KS2) ♦ ( 1 . -PLC ) *  IP  1 I » JK-1.RS2) 

1  -2.*P  (  1 .  JK,Kb2)  )  )  /  ( 1 .  ♦P1_C> 

IF  (J.EO. JR)  GO  TO  184 
p?>xUI*PA*U2*P6*D3*PC 
GO  TO  118b 

18&  P5=01*PA*U2*P8*03*PC 

PAA  =  P  (  I  »  Ji\-1  »KS)  *3,*  (PA-P  (  I  t  JK»KS>  ) 

P0B  =  P ( I f  JR” 1 t  K  S 1 ) *3.* (P8-P ( j . jRtRSl) ) 

PCC  =  P ( I ♦ JR-1 *N5?) *3.* (PC-P( I  * JK»RS2)  ) 

P6=0l *PAA*D2*PBB*D3*PCC 
118S  IF  (I.Nt.lM)  GO  TO  ?30 
P3=0.0 
P4=0.Q 
GO  TO  240 
230  13  =  1 ♦  IP 
1 4  =  I  3* 1 P 

P  3  =  0 1 *P ( I J» J*RS) *D2*P< I3» J.KS1 ) ♦Q3*P < 13* J»KS2) 
P4  =  01*P  (  14.  J.KS)  ♦D2*P<  I4*J»KS1>  ♦l)3*P  ( 14 1  J»KS2) 
240  OENOM=CAPR*A1 ♦CAPL*H1 ♦CAPM*C1 
T 1=-CAPK* ( A?*P3*A3#P4) 

T?=-CAPl* <d?*PS*R3*Pf>) 

T3=-CAPw* (C2*P1*C3*P?) -CAPP 
PbUPF  = ( T 1*12* 13) /DENOM 
P01=F l*P50PF*t2»Pl*F 3*P? 


P01=F  l#P50PF*t2»Pl*f  3< 
P02=3.*pDl-3.*Pl *P? 

IF  (M.EU.2)  GO  TO  250 
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I)L1  (  I* J)=PD1 
0L2<  T ,  J)  =P02 
PWPL  (  I  *  J  >  =PSU*F 
RETURN 

2Bn  D01 ( I ♦ J) =  P!>1 
002 ( i . j) =pd2 

PWBU {  I  .  J) =PSUHE 
RETURN 


PLANFORM-SlOt-fcl>GE  1 M  AGfc  -PO  I  NIT  POIENTIALS 

ENTRY  SURESC 
|ma=IM-1 
OK  =  0E  T  A 
OEL=DELS  ( i ) 

IF  (I.GT.ITE)  GO  TO  300 
Orl  =  -()Xl 
IP  =  -1 
GO  TO  310 
300  OH® ♦OX I 
IP  =  M 

310  FP  =  F  (  I ) 

IF  (I, LT.iLE.OR. I.GT.ITE)  GP  =  GF 

IF  (I .GE.ILE.AND.I.LF.I  IF)  6P=G*T 

GSt T  =  (  (OK-OE'L)  *CAPG  < I  »  J)  ♦DFL*CAPG  (  I  ,  J-l )  )  /OK 

GT A=C APG T 1 ( J«Kw) ♦CAPGT2 ( I  * J) 

GTB  =  CAPGT 1 < J-l «KW) »CAPGT2 ( I  * J-l ) 

GTSET=( (OK-OFL) *GT A ♦DEL*GTF) > /DK 

A  1  =  1 , 5/OH 
A2=-2 • O/OH 
A3=0.5/OH 

Bl=- <0K*2.*0EL) / (DEL* (Ok ♦DEL) ) 

B2= (OK*OEL) / (uEL*OK) 

B3= “DEL/ (OK* (OK ♦DEL ) > 

E1=2.*QK*UK/ ( Of L * (OK ♦DEL ) ) 

E2=-?.* (OK-OEL) /DEL 
E3=<OK-UEL)/(UK*OE'L> 


12=1 1  *  I P 

IF  (I.EU.1M)  11=1 

IF  (I.EU.lMA.OR.I.EO.tM)  !  2= I  1 
J1=J*  l 

P1=P ( I »  J  *  Kw ) 

P2  =  P  ( I »  J  1  -  K* ) 

P3= ( ( OK ♦OtL ) *P ( 1 1 ♦ J ♦ K  w ) — OEL  *P ( 1 1* J1 *KW)) /OK 
P4  = ( <  DK ♦OtL ) *P ( I ?  *  J»K  w ) -DEL*P ( 12* Jl »Kw) ) /OK 
PS  =  (  ( f)K ♦DtL )  *P  (  1  ♦  J.Kw-1 )  -DEL*P  ( I » J1  *Kw-l  )  )  /DK 
P6=( ( OK*OtL ) *P ( I « J » K  w  ♦ 1 ) -DEL*P ( I ♦  u  1 *  K  w  ♦ 1 ) ) /OK 
GG= ( (OK ♦OtL ) *GAMMA ( J) -OEL«GAMMA ( J1 ) ) /OK 
IF  (I.GT.ITE)  P5=PS*GG 

SCl=FP*GStT 
SC2=GP/BSPAN 
SC3=H(KW) *GTSET 
IF  (I.EU.1M)  SC  1 =0 • 

OENOM  =  SC  1  A 1  ♦SC2*H2 
T1=-SC1*  < A2*P3*A3*P4) 

T2=-SC2*(b2*Pl*B3*P2) 

T3=-SC3* (Ph-Pb) / <  2 . *0/FT  A) 

PSURF  = ( T 1 ♦ T2* (3)/DFN0M 

OS  1  (  I  ) =£1*PSURF ♦E2*P1*E3*P? 

0S2 ( I ) =3.*0S1 ( I ) -3 . #P 1 *P2 
RETURN 


WING/rtOOY  SYMMETRY-PLANE  POTENTIALS 


79 


nnonon 


410 


420 


470 


480 


490 

500 


ENTRY  SYMHC 
KwX=KW- 1 
K*P=KW*  1 
KMX=KM- 1 
DO  500  11=2,1* 

GSET=CAPG< I  1*1) 

FC1=GSET*F  (II) /OXI 
FC2  =  G  ( 1 )  /  (RSPAN*IUTM 
KL  =  KH0DL  <  I  I  ,  1  > 

KU=KH0DU (II , D 

IF  (  1 1  .EG*  IM)  GO  TO  4/0 

00  410  K  =  2 » F  R a 

IF  (KL.NF.O.ANO.K.GT.KL)  GO  TO  410 
CAPGT=CAPGTl  (  l  ,K ) *C APGT2 (11*1) 

FC3=CAPGT*n (K) /DZFTA 

IF  (GSET.LE.O.)  PSYM ( 1 1 , K )  =  ( FC 1 *PS YM (1 1 -1 *«> -F C2*P U I  *  1 *K ) 

1  *FC3*PSYM ( 1 1 *K-1 ) )/(FCl-FC2*FC3> 

IF  (GSET,oT,0.)  PSYM  <II*K)s(-FCl *PS  YM ( 1 1  ♦  1  »  K ) ~FQ2*P (II»J»K) 

1  ♦FC3*P5YM(  II.K-l)  )/ (-FC1-FC2*FC3) 

CONTINUE 

00  4  ?0  K  =  M*P,KMX 
KK=KMX-K*KWP 

IF  (KU.NE.O. ANO.KK.LT. KU>  GO  TO  420 
CAPGT=CAPGTl ( 1 ,KK) +C APGT2 ( 11*1) 

FC3  =  CAPGT*H(KK)  /02F.TA 

tF  (GSET.LE.O.  )  PSYM ( l I*KK)  =  (FC1*PSYM(JI»1»KK)-FC2#P ( II » 1 »KK) 

1  -FC3*PSYM< 1 1 «  kK* 1 ) ) / (FC1-FC2-FC3) 

IF  (GSET.GT.O.)  PSYM  (  II  ,KK  )  =  ( -FC  1  *P5YM  (  1 1  ♦  1  »  KK  )  -FC2*P  (I  I  »  1  »I^K ) 

1  “FC3*PSYM( 1 1  * KK *  1 ) ) / (-FC1-FC2-FC3) 

CONTINUE 

IF  ( 1 1 *6E. IN05E)  GO  TO  500 
CAPGT  =  CAPGT 1 < 1 » F.W ) *CAPGT2 ( 1 1 » 1 ) 

FC3=0.5*CAPGT*H(KW) /DZETA 

PSYM(II,K*)  =  <FCl*PSYM(n-l»KW>-FC2*P<II.l»KW)-FC3MPSYM(II,KW*l> 

1  -PSYMUI.KW-1  )  )  )/ (FC1-FC2) 

GO  TO  500 
00  480  K=2,KL 

CAPGT sCAPGT 1 ( 1 ,K) ♦CAP6T2 (1 1  *  1 ) 

FC3=CAPGT*H(K) /OZETA 

PSYM ( IM,K) s(-FC2*P < IM, 1,K) *F C3*PSYM ( IM,K-1 ) ) / (-FC2*FC3) 

00  490  K=KU,KMx 
KK=KMX-K*KU 

CAPGT  =C APGT 1 ( 1 , K  K ) *C  APGT  2(11*1) 

FC3=CAPGT»H(Ka) /DZETA 

PSYM (IM,KK)=(-FC?*P(IM,1 ,KK)-FC3*PSYM (  JM*KIW  ) ) / (-FC2-FC3) 

CONTINUE 

return 

ENO 


SUBROUTINE.  TRIC0E(K1,*2,I,  J) 

CALCULATES  FINITE-DIFFERENCE  COEFFICIENTS  AT  GRID  POINTS  ALONG 
The  SEGMENT  K1  TO  K2  OF  LINE  (I»U>. 


MMON/CA 


COMMON/COEF  / 


C APH  (49*19)  *CAPHT1  (19,25)  ,CAPHT2  <49,  19)  ,CAPI8U9>  ♦ 
CAPQB (49,19) ,CAP IT (49,19) 

A1 (49) ,A?(49) ,A3(49) ,A5(49) ,81(19) *B2(19) ,B3(19) , 
B5( 19) ,B7 (19) ,CI (25) * C2 ( 25 ) *  C3 ( 25) ♦ C5 ( 25) , C 7  < 25) , 
D1,D2,E1,E2,FJ,F2 


80 


nnnno 


10 


20 


30 


COMMON/CONST/  ALPhA.ZmaCH.DPLIM.EPSI .WE.KG.BSPAN 
SNA.TDETA.TDXI .TOZETa.DETA.DXI.OZE 
ITAIL*ITE*  jm»  jroot .JTIP.KM.Ka.JBl * 
PLCH2.PLCB3 

COMMON/COUHO/  CHOROU9) .DCDY(19> .DTDY<19) .DXLEDY 
.FLE.FN.FTE.G( 19) * GAMLE ♦ GAMN ♦ GAMF , 
T AU (19)  .  X ( 49 *  1 9 ) ♦ XCAP (49) «XLE (19)* 
2 CAP (25) • ZLE (19) .UZlEDY(19> » ALPHA T 

common/solvo/ 

JMAXI ♦ JMAXT.KMAXl»KMAXT.NSUP.P<49. 


.AI2.CSA.RX1.RX 2. 
TA.ILE.IM. I NOSE* 
JB2. JB3.PLCB1 . 


PNEW (25) ,PN0SE»PSYH(49.25) »PT1  (19. 
PT3 (19*25) .PTE  1 (19) »PTE2(19) .P*8L( 


(19) *ET  A ( 19) » F ( 49 ) 
GAMTE.GAMmT  *H(2&)  * 
Y ( 19) ,ZETA (25) *GF, 
(19) » OAOY ( 19 ) »GWT 

tSliJVHISf}”* 

19.25) .PLE1 (19) ♦ 
25) .PT2U9.25) * 
49.19) *PW0U (49* 19) 


CHECK  P.O.E.  TYPE 


FP=F (I) 

GP*G ( J ) 
CJ»CHORO ( J) 
TC*TAU(J)#CJ 
CBG*C APG ( 1 ♦ J) 
CU1*FP/CJ 
CU2-=CAPGB  ( I  .  J) 
:v|*C8G*FP 


CV2=GP/8SPAN 


^WlafP/aJ-CJ) 


_.»FP*CbG«CdG 
CAP3=-2.*FP*CdG/CJ 
C8PUGP/  (BSPAN*8SPAN) 

CCP 1 *CU2*CU2 

CCP2=1 ./ ( TC*TC) 

CCP3*-2.*CU2 

CCP4a-2./TC 

CCP5=-2.*CU2/TC 

CDP1=2.*FP*GP*CBG/BSPAN 

CDP2=-2.«FP*GP/ (CJ«8SPAN) 

CEP1»2.*GP/BSPAN 

CEP2»-2.*GP*CU2/BSPAN 

CEP3=-2.*GP/ ( TC*BSPAN) 

CFP1«2.*FP*CU2/CJ 

CFP2a2.*FP*CBG 

CFP3=-2.*FP/CJ 

CFP4a-2.*FP»CU2*CBG 

CFP5*-2.*FP*CdG/TC 

CFP6*-2.*FP/ <  TC*CJ) 

‘  *  - ( J)  *CJ 


t JP 1 *2 . *CAP 
'JP2=-CAPH( 

3P3=«2.*CAP 


AP 


J|jP4*2i* 
CJP5a-TC 
CJP6a-CU2*TC 


♦J)*Cv 

T  1 1*  J  j 


8  ( J)  * 


*TC 

C 


JKaJM 
JKP= JK ♦ 1 

IF  ( J.EU.JB1 .OR. J.EO. JB1M ) 

IF  (J.EQ.JB2.0R.J.FQ.JB2.1) 

IF  ( J.EO. JB3.0R.J.EQ.JB3.1 ) 

GO  TO  90 

JKsJBl 

PLC=PLC01 

GO  TO  80 

JKeJB2 

PLCaPLCB2 

GO  TO  80 

JK»JB3 

PLC»PLCB3 


GO  TO 
GO  TO 
GO  10 


10 

20 

30 
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noon 


an  jkf= jk *  1 
90  CONTINUE 

00  700  K=K1,K2 
HP*H (K) 

GTP=CAPGTI(J.K)  ♦CaPGT<?(  I.J) 

HTPsCAPHTI ( J,K) ♦CAPHT2 ( I .  J) 

P1  =  P< I»1.J*K> 

P2«P ( I • J«K) 

P3*PT2(J»M 

P5=Pn»J»K-)l 

P8*Pil« 

IF  (J.EQ.U  P 1 1  *PS  YM  ( I  *  K ) 

IF  (J.NE.D  PH»PT1 
P16*P(I*JMtK> 

If  (J.EQ.JK)  P16=<2.*P16M1.«PLC)*<P11-2.*P2)  >/U.*PLC> 

IF  ( J«EQ»UKP)  P11=(PLC*P16-2.*(PLC*P2-Pil) )/(2.-PLCI 
UP=CSA*CU1*<P1-P3) /T0XI*CU2*HP* l  PB-P5 ) /TOZET A 

VP=CV1* (P1-P3) /T0XI*CV2* (P16-PI1) /TOETA»GTP*HP* (P8-P5) /TDZETA 
WP*SNA*CWl*HP*  <Pft-P5) /TDZETA 

{£  (i^kf  I*W^AnS?kInE^*1>  GO  TO  105 
UP=UP-CU2*HP*GAMMA ( J) /TDZETA 
VP*YP-GTP*HP»GAMMA  < J)  /TDZETA 
WP*WP-C*1*HP*GAMMA ( J) /TDZETA 
105  UPSQ*UP*UP 
VPSQsVP*VP 
WPSU*WP*WP 
QPSQ«UPSQ*  VPSd* WPSQ 
TQ*A 12*0. 20-1 ,20*QPSQ 

apsg*to*qpsq 

IF  (TQ.LT.O.)  GO  TO  200 

DEFINE  POINTS  OF  ELLIPTIC  (SUBSONIC)  COMPUTATION  MOLECULE  AND 
COMPUTE  TKI-OI AGONAL  COEFFICIENTS 

P3NsP ( 1-1 ♦ J»K) 

IF*U  IeWtE.ANO.KI.EO.KW*!)  P4*P4* GAMMA  ( J) 

P6N=P(I-1»J,K-1) 

P7*P(I*1»J*K*1) 
p9N=p<  i-i  »j,k*d 
P15«P(l*I*J*l*K) 

Pl7NaP(I-l,J*itK> 

P18*P(ItJ*l»K-l) 

P19=P(I.J*l.K*l) 

PI 1N*P ( I  * J-l «K) 

P12N*P(I-I» J-I.K) 

P13N»P(I.J-1 ,K-1) 

P14N=P<I.J-1»K*1) 

!f  ( J»EQ*UKi  GO  TO  112 

o  mibJKP)  Go  To  116 

1 0*PS VM ( I ♦ 1 1 K ) 

PllN*PSYM<I,K) 

Pl2NaPSYM<I-l*K) 

pi5n»psym(i,k-i ) 

PI4NePSYM(f*K^t) 

GO  TO 
P15* (2 
P'7N*I.  t(1. 

p  .e«(2»*pi8«(ia-  ...  ..  ....  _ 

p  9*(2,»P19* (I.-PLC)*<P14N-2,*P8) )/ 

GO  TO  120 

pI0«(PLC*P15-2.*(PLC*P1-P10) )/ (2.-PLC) 


no 


1 12  P 


z'.Ux*  ♦(1»-PLC)*(P1Q-2»#P1)1/(1 •♦PLC) 
12.*Pl7N*a.-PLC)*tP12N-2,*P3NI  )/(I,*PLC) 
2.*P18*  1.-PlC)*<P13N-2,*P5) )/(1.*PlC) 


116 


120 


H:A8 


P11N*(PLC»P16-2.*(PLC»P2-P11N> )/(2.-PLC> 
P12N*(PLC#P17N-h,* (PlC*P3N-P12N) ) / (2. -PL 


PLC) 


pl3N»<PLC*Pl8-2,*(PLC*P5-P13(ii)  )/(2.-PLC) 
P14N=(PLC«P19-2,*<PLC*P8-P14N) >/<2.-P LC> 

TU*APSO-UPSO 

TVbAPSQ-VPSQ 


82 


nnnn 


TWaAPSO-wPSO 

CAP*CAP1*TU*CAP2*TV*CAP3«UP*VP 

£cPs^Cc4T*TU*GTP*<TV*GTP*CCP3*UP*VP*CCP4*VP*WP) *CCP2*TW 
1  ♦CCP5*UP**P) *hP 

C0P*C0P1 * TV*C0P2*UP*VP 

CEP=  <CEP1*TV*GTP*CEP2*UP*VP*CEP3*VP*WP) *HP 
CFP= (CFP1*T U*CFP?*TV*GTP *UP* VP* <CFP3*GTP*CFP4> 

1  *WP* <LFP5*VP*CFP6*UP>) *HP 

CJP= (CJPI*UP*VP»CJP2*TV) « <UP-CSA*CJP6*(WP-SNA) ) 

1  ♦ (CJPi*UP*VP*CJP4*VP*WP*CJP5*rv*HTP)*(WP«SNA> 

ACAP(K) *CCP*C3(K) 

8CAP(K) *CAP*A2< I ) *RX1 ♦CSP*B2 ( J) *CCP*C2  <K) 

CCAP(K>«CCP*C1 (K) 

OCAP(K) *CJP-CAP*(A1 (I ) *P1*A2( I) *RX2*P2* A3 ( I ) *P3N) -CBP* 

1  <B1  <  J)  *P  16*83  ( J)  *P1  IN)  -CDP*(Dl*(P15*P12N)*D2*(P10*P17N)  )  - 

2  CEP*  <£1*<P19*P13N) *E2* (P18*P14N) ) -CFP* (F 1* (P7*P6N> *F2* (P4*P9N> ) 
IF  (K.EQ.Kl)  UCAP (K) eDCAP (K) -ACAP (K) *P5 

IF  (K.EQ.K2)  DCAP (K) «DCAP (K) -CCAP (K) *P8 
IF  U.LE.ITE.OR.J.LT.JROOT)  GO  TO  700 

IF  (K.EQ.KM)  0CAP(K)*DCAP(K)*CCP*CI<K>*GAMMA<J)-CEP*E2* 

1  (Gamma  (j*  l  > -Gamma  (j-m 

IF  (K.EQ.Kw*1{  UCAP(K)«0CAP<K)-CCP*C3<K)*GAMMA(J)-CEP*E2* 

1  (GAMMA(JM)-UAMMA(J-1>  > 

GO  TO  700 

DEFINE  POINTS  OF  HYPERBOLIC  (SUPERSONIC)  COMPUTATION  MOLECULE  AND 
COMPUTE  TRI-DIAGONAL  COEFFICIENTS 


200 


NSUPsNS UP* 1 
P3N*P ( 1-1 ♦ JfK) 

P4*P ( I ♦  I  • U  * K-  1  ) 

IF  (I.EJ.1TE.ANO.K1.EQ.KW*!) 
P6NaP(I-l«JtK-I) 
P7*P(I*I»U.K*1) 

P9N=P ( I-l . J*K*1  ) 
P15«P(I*1*J*1«K) 


P4*P4« GAMMA ( J) 


t-i# J*1 »K) 
(  .J*  .K-  ) 

(!*j*1«k*i) 

_  «EU«  1 )  GO 

lIBIWfi 


P 1 7NsP (I-l 
P I 8aP ( I »  J* I  *  K" 

P I 9*P (  I  »  J  ♦  I  ♦  K  < 

£F  (J.EU.D  GO  TO 


270 


P 1 3N*P (IfJ— 1»K  —  1) 
P14N*P(I» J-1»K*1) 


270 


272 


276 


280 


* 

i 


IF 

F 

0 

P 

P 

P 

P 


(J.EO.UK)  GO  TO  272 
(J.EU.JKP)  GO  TO  276 
TO  280 


0  =  PSYM ( l *  1 tft) 

1N«PSYM ( I »K ) 

2N=PSYM( I-l  * K ) 

3NaPSYM ( I *K— 1 ) 

P14N*PSYM ( I »  K ♦ 1 > 

GO  TO  280 

P15»(2.*Pi5* U.-PLC)»(P10-2.*Pl) )/(l.*PLC) 
P17N*(2.*P17N* (1 .-PLC) * (PI2N-2,*P3N) ) /(l.*PLC> 
PI8*(2.*P18* ( I,-PLC)*(P13N-2.*Pb> )/( 1 ,*PLC) 
g49T02*8P19* ( l*“pLC> *(P14N-2.*P8) ) /< l.*PLC> 

?l|M?Se8»fi5s:i7b£e5»5»?U«f(»Kc» 

P12Na(PLC#Pl7N-2.*(PLC*P3N-Pl2N) ) / (2.-PLC) 
PI3N«(PLC*P18-2.*(PLC*P5-P13N) )/ (2,-PLC) 

P14N* (PLC#P19-2.* (PLC*P8-P14N) >/ (2.-PLC) 
CONTINUE 


P20*PT3 ( J  *  K ) 
IF  (K.EQ.M) 
IF  (K.NE.KJ 
IF  (K.EQ, 

IF  (K  »NE i 


i  • r'  *  i 

i:S|l 

,.K2) 


P21*P5 

P21=P(I«J»K-2) 

P22»P8 

P22aP(I»JfK*2) 
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IF  (J.EQ.l)  P23*2.*Pll-P2 
IF  (J.EU.2)  P23*PSYM( I ,K) 

(j.GE.i)  r 


IF 


P23*PT1 (J-2»K) 


CPP  = 
CQP* 
CRP  = 

*CSP* 

CTP* 

CVP» 

i 

cwp* 

LUCP2 

UCQ2 

UCR2 

UCS2 

UCT2 

UCV2 

UC*2 

FEQ 

CLP 

CMPs 

CNPa 


CAP1*UPSQ*CAP2*VPSQ-CAP3*UP#VP 

CBP1*VPSQ 

(CCP1#UPSJ*GTP*  <  VPSQ*GTP-CCP3*UP»VP-CCP4*VP*tfP>  ♦CCP2*#PSQ 

-CCP5*UP**!P>  *HP 

CDPI»VPSQ-CDP2*UP*VP 

(CEPl*VPSQ*GTP-CEP2*UP*VP-CtP3*VP*WP) *HP 
<CFPI*UPSQ*CFP2*VPSQ*GTP-UP*VP*  <CFP3*GTP*CFP4> 

-*P«  (CfP5»VP*CFP6*UP)  )  *HP 

(CJP1#LIP*VP-CJP?#VPSQ>*  <UP-CSA-*CJP6*(KP-SNA)  ) 

♦  (CJP3*UP*VP>CJP4*VP*WP«C  JPS*'/PSQ*HTP)  *  (WP-SNA) 

»QPSQ* (CAP  1 ♦CAPE )-CPP 
»OPSQ*CBPi-CQP 

»QPSQ*(CCP1>GTP«GTP4CCP2>  *hp-crp 

*QPSQ*CDP1-CSP 

aQPSQ*CEP 1*6TP*HP-CTP 

sQPSO* (CFP1*CFP2»GTP) »HP-CVH  .  „  „  v  „  „ 

*QPSQ# < -C JP2* (UP-CSA*CJP6* I WP-SNA) )-CJP5*HTP*(WP-SNA) ) -CWP 
FP*UP/ (CJ*DXI) 

FEQ*  <  UP/CQ ♦ VP*CB6) 

FEQ* VP/BSPAN 

FEQ* (QP*CU2*GTP*VP*WP*CW1 ) 


1 

2 


TQ*1.- 

ACAPT* 

BCAPT 


CC  APT 
OCAPT* 


1 

§ 

4 

1 

1 


OCAPT  a 


QPSU/APSQ 
UCP<?*C 3  ( K  ) 

;m:i  lumnnmsfffi  wnuau*' 1  ’ 

-EPSI* (CLP*A3(I)*DXI ♦CMP*B3 ( J) *DETAJ 
UCR2*C1 (K) 

-UCP2*<A1 <I)*(P1-P2) ♦A3(I)*P3N) 

-UCQ2*(B1  (J)*(P16-P2)  ♦B3(Q>*PUN> 
-UC52*(01*(PI5^P12N)  *02*  EI2*£PN) 
-UCT2*(E1*(P19»P13N>  *E2* (P18*PI4N> > 

-UCV2* (FI* (P7«P6NI ♦F2*(P4*P9N> )-UCW2 
DCAPT-TQ* 

(CPP* (-A3 ( I ) * (P2*P3N> -A5  < I >  * (P3N-P20)) 
♦CQP*(-B3(J)*(P2^P11N)-BSU)*(P1IN-P23)  ) 

♦  4,*CSP*<01*P12N402*(P3N*PUMI  >  +C»P> 
OCAPT=OCAPT*EPS I  * 

l  ( CLP* A 3 ( I) *0X I* (-P2»P3N*P3>  v 

\  ♦CMP*B3(J) *DETA* (-P2-PI1N*P11) > 

IF  (WP.EQ.O.)  GO  TO  320 
IF  (WP.LT.O.)  GO  TO  330 

ACAP  (K)  a  AC  APT  *TQ*  l-CRP*  (C3(Kl  *C5  (*)  )  *4.*  <CTP*E2*CVP*F2)  ) 
♦£PSI*CNP*C3 IX > *02ETA  _ 

8CAP(K) =BCAPT*TQ*2,* (CRP*C3(K) *2.* < CTP*E 1 *CVP*F 1 > ) 

-EPSI *CNP*C3 (K ) *02ETA 

! aO^APT-TQ*(CRP*<-C3(K) *P2*C5(K) *P2i) 
♦4,*CTP*(E1*P13N+E2*P11N> 

♦4,*CVP*(F1*P6N^F2*P3N) ) 

♦EPSI*CNP*C3(K) *02ETA*(-P2*P5> 

GO  TO  400 


AP(Ki 
AP  ( K  1 


320  ACAP ( K 

^BCAP ( K 
CC AP ( K 

^  DCAP ( K 
1 
2 

GO  TO 

330  ACAP  ( K 
8CAP (K 

1 


»*ACAPT*TQ*(CRP*C3(K> ♦2,*(CTP*E2*CVP*F2) ) 

♦  EPS  I *CNP* HP/ TOZET  A 

)*8CAPT-TQ*CRP*(Ci(K}^C3<«) >  „ 

)*CCAPT*TQ*(CRP*Ci(K) ♦2«*(CIP*EI*CVP*FI) ) 

-tPSI*CNP*HP/lDZETA 
) aDCAPT-TQ* (2.*CTP* (E 1 *P l 3N^E2*P 14N) 
♦2,*CVP*(F2*P6N>F2*P9N) > 

♦EPSI*HP* (-PB^PS) /TOZETA 
400 

) sACAPT 

) sBCAPT*TQ*2,*(CRP*Cl (K) *2, •  ICTP*E2*CVP*F2) ) 

♦  EPSI*CNP*C.l  <*)  *DZETA 
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CCAP  (K) sCCAPT ♦ TO* < -CRP* ( C 1 < K >  *C  M K )) *4 .  *<  CTP<*E  1  ♦CVP*F  1 ) ) 
1  -£PSI*CNP*C1 (K) #DZETA 

OCAP  (K)  cOCAP  T-TQ*  (CPP*  ( -Cl  (K)*P2*C7(K)*P22) 

1  ♦4.*crP*<E2*PUlM*El*PllN> 

2  ♦4,«CVP* <F2#P9N*FI«P3N) > 

3  ♦£PSI#CNP»C1 (<)*0ZETA* (P2-P8) 

400  CONTINUE 

IF  (K.EQ.K1)  0CaP(K)bDCAP(K)-ACAP(K)*P5 
IF  (K.EU.N2)  UCAP(K)»0CAP(K)-CCAP(K)»P8 
700  CONTINUE 
RETURN 


750 


COMPUTE  TRI-OIAGONAL  COEFFlCIENrS  IN  THE  THEFFTZ  PLANE 

ENTRY  TRIT 
GP=G( J) 

TC=TAU<U)#CH0R0(J) 

CV2»GP/8SPAN 

C*1=1./TC 

00  900  K*K1,K2 
HP*H (K) 

GTPaCAPGTl (Ut*) ♦CAPGT2 ( IM, J) 

HTP*CAPHT1 (J«K) *C APMT  2 ( I M  t  J  > 

P5»P(IM,JtK-l) 

P8«P( IM,J»K*1) 

IF  (J.EQ.l)  P1I=PSYM( IM*K) 

IF  (J.NE.l)  PI1=PT1 (J-I.K) 

P16»P(IM.UM*K> 

VP»CV2#(P16-PU  )  /TDETA-»GTP*HP#  (P8-Pb)/T0ZETA 
WP»SNA*CWl#HP* (P8-P5) /TOZETA 
IF  (K.NE.^R.ANO.K.NE.KW^l )  GO  TO  750 
VP»YP-GTP*HP*GAMMA ( J) /TOZETA 
WP=RP-CW1*HP»GAMMA ( J) /TOZETA 
CONTINUE  ^ 

TV=AI2-VP#VP 
CBT3TV*GP/ (BSPaN*8SPAN) 

rrTaHP* ( (AI2-RP*WP) / ( Tr*Tc) ♦  GTP* ( Tv*gTP-2.*VP*WP/Tc> > 
CET»2. #GP#HP# (TV*GTP-VP*WP/TC)/B5PAN 
CJT*TC* (WP-SNA) * (2,*VP*WP*CAPIB( J) -TV«HTP) 


•u 


800 


850 


P2*P< IM.J.K) 

P18»P ( IM» J*1 «K-1 ) 

P 1  9*P ( I M  t  J* 1 **♦ 

IF  <J.EQ.D  GO 
P 1 1 NsP ( I M  *  J- 1 »  K ) 
P13NoP( IM»U-1 .K-l  ) 
P14N«P ( IMi J-l *K*1 ) 
GO  Tn  850 
PUN»PSYM<  IM»K) 
P13N*P$YM IIMiK-1 ) 
PI4N=PSYM ) 


800 


ACAP (K) *CCT*CJ (K ) 

BCAP(K) *CdT*82( J)*RX1*CCT*C2(K) 


ss 


1 


AP  (K) 


IF 


_UT-C8T* ( B 1 ( J) *P16*B2(U)*RX2*P2*B3(U)*PI1N) 
-CET* (El* (P19*P13N) *E2# (P18*P14N) ) 

(K.EQ.K1)  OCAP(K) =DCAP(K)-ACAP(K) *P5 

_  -  -  -  ••  ‘  “  " 


1 


(K.EQ.K2)  0CAP(K>=0CAP(K>-CCAP(K)*P8 
(J.LT.UROOT)  GO  TO  900 

. .  .  . "  ~ "  T  *  C 1  (K)*GAMMA(J) 


900 


1  -CET*E2# (Gamma ( j* i ) -gamma  < J-l ) ) 
CONTINUE 
RETURN 
END 


CT*C3(K)#GAMMA(J) 
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SUBROUTINE  INVERT (K1 «K2) 

INVERTS  TRI-DIAGONAL  MATRIX  (VARGA,  MATRIX  ITERATIVE  ANAL..  P.I95) 


DIMENSION  GG<25) ,WW(25) 


COMMON/ABCO  / 
COMMON/SOLVO/ 


AC AP ( 25  ) ,BCAP<25) ,CCAP<25> ,0CAP<25) 

DL1 (49» 19) ,0L2<49» 1 V) .DSl  (49) .052(49) »DU1 (49.19) ♦ 
DU2 (49. 19) . DPMAX . UPMAXT  » GAMMA ( 1 9 ) » I  MAX  I . I TERT » 
JMAXI.JMAXT.KMAXI.KMAXT ,NSUP«P (49.19.25) ,PLE1 (19) ♦ 
PNEW ( 25) .PN0SE,PSyM<49,25> ,PT 1(19.25) .PT2( 19.25) » 
PT3( 19,25) , PTE  1 ( 19) ,PTE2(19) »PW8L <49* 19) »PW8U(49*1 


CALCULATE  COEFFICIENTS  OF  REDUCED  SYSTEM  OF  EQUATIONS 

WW <*1  ) =CCAP (Kl ) /RCAP (Kl) 

G6  <  K 1 ) *UCAP (Ki ) /BCAP (KI) 

K  1  P=K 1*1 
00  100  K=N1P,*2 
KZ= K-l 

OtNOMsBCAM (K)-ACAP(K)*Wd<K7) 

WW <K ) =C CAP (K ) /DENOM 

GG (K ) = (OCAP (K ) -ACAP (K ) »GG (K^) ) /DENOM 

INVERT  REDUCED  MATRIX  BY  HACK  SUBSTITUTION 

PNEw (K2) =uG (K2) 

KI)  =  K2-Kl 

00  200  K=l, KD 

KZ=K2-K 

PnER(KZ)=UG(KZ)-WW(KZ)*PNEW(K7*1) 

return 

END 


SUBROUTINE  MAXI (Kl ,K2, I , J) 


determines  location  and  absolute  value  of  the  maximum  change  IN 
potential  between  the  last  two  computed  iterates. 


COMMON/SOLVO/ 


DL1  (49,19)  ,DL2  (49, 19)  ,DSl  (49)  ,DS2<49)  ,DU1<49,19)  ♦ 
DU2 ( 49. 19) ,0PMAX,DPMAXT,GAMMA(19) , IMAXI ♦ I TERT , 

JMAXJ ,JMAXT,KMAXi»KMAXT,NSUP*P(49.19.25) ,PLE1 (19) ♦ 
PNEW (25) ,PN0SE,PSYM(49,25) »PTi (19,25) *PT2( 19.25) ♦ 
PT3( 19.25) ,PTE1 (19) ,PTE2(19) .PwBL(49.19) ,PW8lM49*l' 


00  100  KsKl,K2 
OELP=PNEW(K)-P<I,j.k) 
OP=ARS (OELP) 

IF  (DP.LT .OPMAX)  GO  TO 

OPMAXsOP 

I  M  A  X  I  s  I 

JMAXTsJ 

KMAXIsK 

CONTINUE 

RETURN 

END 
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nr>c->  nnficnr, 


SUBROUTINE  ARRAYS(NPRINT) 

OUTPUTS  SOLUTION  ARRAYS  EVERY  NPRINT  ITERATIONS.  IF  NPRJNT.EQ.Q, 
ONLY  THE  CIRCULATION  IS  WRITTEN. 


COMMON/CONST/  ALPHA»2HACH»DPLIM*|RSI»wE.we*BSPAN,Al2*CSA.RXI.RX2. 

*  SNA.TOETA.TDXI «T02ET  A»OETA»OXI .OZETAt ILE » IM, I  NOSE. 

*  I  TAIL  « ITE, JM, JRQOT.UT IP ,KM» KW, US l « J82* J83.PLC8 1 ♦ 
pirQPiPic83 

COMMON/SOLVO/  DL1  ( 49 , 1 9 ) , DL2 ( 49 » 19) »QS1 (49)  t DS2 < 49) »OUl (49* 19) ♦ 

*  OU2  U9. 19)  ,OPMAX.OPMAXT,GAMHA ( 19) ♦ I MAX  I « ITERT. 

*  JM  AX  I  .  JMAXT  ,KMAXI »*MAXT  »NSUP  »P  $  49 . 19  ♦  25)  ,  PLE I  ( 19 )  » 

*  PNEW (25) .PNOSE  »PSTM< 49,25) *PTl(I9»25)jPT2(l9»25)* 

*  PT3( 19.25)  .PTE  1 { 19) .PTE2 (19) »PWBL (49, 1 9) ,PW8U<49» 19) 


L 


C 


C 


c 


c 


WRITE  (6.161) 

WRITE  (6.171) 

WRITE  ( b .  1  7 1 ) 

WRITE  ( 6. 1  ft 1 ) 

WRITE  (6.161) 

WRITE  (6,131)  (GAMMA (J) ,J=1 »JM) 

IF  (NPRINT.EO.O)  GO  TO  225 
WRITE  (6,161) 

WRITE  (6,141) 

UO  ISO  J=i.JM 
WRITF  (6,142)  J 
UO  145  1  =  1.  IM 

145  WRITE  (6,146)  I , <P< I ,J,K) ,K=1 .KM) 
150  CONTINUE 

WRITE  (6,161) 


105 

110 


WRI  IE 
UO  110 
WRITE 
00  105 
WRITE 


1 


(6.101) 
J=1 . JM 
(6,102) 
1  =  1  ,  IM 
(6,106) 


CUNT INUE 
WRITE  (6,161) 


J 

I , PWBL l I ,US ,0L1 ( 1 »U>  »DL2 ( I . J) *DU2 1 1 , J) ,UU1 ( I, J) , 
PwRU ( 1 , J) 


WRITE  (6.121) 

UO  125  1*1. IM 

125  WHITE  (6,126)  I  ,r>Sl  ( I )  ,DS2  U  ) 
WRITE  (6.161) 


?20 

225 


WRITE 
UO  220 
WRITE 
WRI  TE 
WRITE 
WRI  TE 
WRI  If 
RETURN 


(6,191) 
1  =  l  ,  I M 
(6,14b) 
<6. Ibl ) 
(6,171) 
(6,171) 
(6,161) 


I  ,  <P5YM(I , K )  , K  s  J ,6M) 


101  FORMAT 
1 

1  02  format 
1 

106  FORMAT 

121  FORMAT 
1 

126  FORMAT 

n\  format 
i 

141  FORMAT 
1 


(  1 H-/// l x 40HDUMM Y-P0 1  NT  POTENTIALS  ON  VERTICAL  LINES/ 
1X2SH(PLANE-HY-PLANE  SPANWISE)) 

(  /  // 1 X 1 OHPL ANE  J  =  I3//3XlHl,l6X4HPWBL,12X3hOLl,12X3HOL2, 
1  7X3HUU2, 12X3H0U1 .1 1 X4HPW6U/ ) 

( I4»2<5X, 3E 15.4) ) 

(1H*///1X3BHOUMMY-POINT  POTENTIALS  ALONG  SIDE  EDGE// 
3XIhI,17X3HDS1,1?X3H0S2/> 

( I4,5X,2£IS%4) 

( 1H-///1X34HCIRCULATI0N  DlSTRlBUT ION,  GAMMA(U)// 

<  4X  ,  l  JE I  f) ,  3)  /// ) 

(  1H-///1X?6HP( 1 , J,K>  :  POTENTIAL  ARRAY/ 
IX43H(PLANE-RY-PLANE  SPANWISE,  I  DOWN,  6  ACROSS)) 
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non  o  nonnn 


142  F ORM  A  T 
146  FORMAT 
161  FORMAT 

171  format 
1 
2 

181  FORMAT 
191  FORMAT 


(///1X10HPLANE  J  =  13/) 

< 14/ (4A, ]  JEJ o. J)  ) 
(//////) 


(  1H-///1X15HS0LUT10N  ARRAYS) 

Ufi-///1X46HPSYM(  I  ,K)  S  IMAGE  POTENTIALS  AT  SYMMETRY  PLANE/ 
1  X 1 8H ( I  DOWN.  K  ACROSS)/) 


SUBROUTINE  CPCOMP 

COMPUTES  THE  WING/BODY  SURFACE  PRESSURE  DISTRIBUTION. 


DIMEN 

DIMEN 


SION  CPL<49, 
SION  ZXLU9. 


19)  *CPUU9»  19) 
1 9 ) » ZXU ! 49  » 1 9 ) 


COMMON/CAPS  / 
COMMON/CONST/ 

COMMON/COORD/ 

COMMON/SOLVO/ 

COMMON/SURF  / 


XLpSSJ|H;iM56SfJiji4Ui!i.«e.BSPAN.Al2,CSA.RXl.RX2t 

SNA.TDETA, TDXI.TOZETA.DETA.OXI.DZETA.IlE.IM.INOSE. 
ITAIL*ITE« JM.JROOT  » JT IP »KM . KW * JB1 * JB2 ♦ JB3.PLCB1 ♦ 
PLCB2 f PLCB3 

uLl (49. 19) »DL2 (49.19) ,DSl (49) .DS2 (49> .DU1 (49. 19) ♦ 
DU2 (49*19) .DPMAX.DPMAXT.§AMMA(19) ♦ IMAX I . I TERT . 

UMAX  I .JMAXT .KMAXl .KMAXT . N5UP  »P(49» 19*25) »PLE1 (19) . 
PNEw ( 25) , PNOSE . PS YM (49.25) .PT1 \ 19.25) .PT2 (19.25) . 
PT3( 19.25) .PTEl (19) ,PTE2(19) .PWBL (49. 19) .PWBU (49. 1« 


KBODU  <49. 19) .ZL (49. 19) »ZU (49. 19) 


♦0.20/AI2)/1.20)**3.5-1.)/0.70 


CPS  T  AR  =  A  1 2* ( ( ( 1 
DO  100  1*1. IM 
DO  100  J=1.JM 
ZXU ( I . J) =0. 

ZXL  (  I . J) =0. 

CPU(I.J)=0.0 
100  CPL ( I . J ) *0 . 0 

OK  =  OE  T  A 
imx*im-i 

00  200  MSURF*  1.2 

DO  200  IsINOSE.IMX 

JE= JPOD ( I ) -1 

00  200  Jsl.JE 

IF  (MSURF. EU. 2)  GO  TO  110 

KP  =  -1 

OL=-DZET A 

UEL  =  -0ELL  I  I . J) 

IZ=IZL (J) , 

ks=krodl ( 1 *U) 

^  ^  ^  ^  ^  ^  ^  | 

CAPGTB* (OL-DEL) * (CAPGT1 ( U. KS ) -C APoT 1 (J.KS*1) )/DL*CAPGTl (J.KS»1) 
CAPGTG=CAPGTB*CAPGT2( I . J) 

GO  TO  120 
110  KP=M 
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c 


c 


c 


■>L=07ETA 
OEL=OELU  <  I  .  J) 

IZ=17U<J) 

KS=NHODU ( 1 ♦ J) 

H  *  B  =  H WBU  (  1  *  J ) 

CAP6T«= (DL-DEL )  *<CaPOt  1  ( J.KS) -CAPGTl  (J.KS-l)  )/DL*CAP6Tl  (J»Kb-l) 
CAPpTB=C APbTH*CAPGT? ( I » J) 

120  IF  (J.Gt.JHOOT)  GO  TO  140 
OH=-OX I 
IP=-1 
GO  TO  180 

140  IF  <I.GE.i2>  GO  TO  160 
ON=-OXl 
IP  =  -1 
GO  TO  180 

160  OH=*OXl 
IP  =  M 


180  A  1 =- 1 ,6/DM 
A2=2.0/0h 
A3=-0.5/0n 
61=-1 .5/0* 

82=7.0/0* 

B3=-0 .5/0* 

Cl=-  ( DL ♦ 2 •  *DEL  )  / (PEL* <DL*0FL) ) 

C2  =  <0L*DEL)  /  <uEL*Ol  ) 
C3=-nEL/(OL*<OL*D£L) ) 

01  =  0.5*  (DL*  DEL)  *  ( 2 .  *OL*0El )  /  ( DL*t)L  ) 
02=-PEL* (2.*DL*0EL)  /  ( l)L*0L  ) 
Q3=0.5*0El* <0L*DEL> /<UL*0L> 


FP  =  E  (  I  ) 

GP  =  6< J) 

CJ=CH0HD ( J) 
CBG=CAPG< I  * J) 

CU 1 =  FP/ C J 
CU2=CAPGR (  I  .  J) 
CV1=CRG*FP 
CV2=GP/8SPAN 
C  W 1  =  1 .  / (TAU(J) *CJ) 


181 


i  18? 


18? 

163 


P2=PW4L ( I « J) 
PZ=P*BU ( I , J) 


I3=I*IP 
I4=I3*IP 
JS=J* 1 
.  J  6  =  J  ♦  ? 

KS1=KS*KP 
Kb2*KSl*KP 
IF  (MSURF.EU.l) 

IF  (MSUPF.EQ.2) 

P 1 =P  < I ♦ J»*S) 

P?  =  P ( I  ♦  J.*S1 ) 

P3aD] *p  l  Ii* Jt*S) ♦D?*p ( 13* JtNSl ) *03*8 ( I3t J»«S2) 

P4  =  I)1*P  (I  4,  J»K5)  *D?*P  <  I  4»  J*KSl )  ♦03*P  (14  t  J»K52) 

Pb  =  Dl*P  <liJ5.*S)*02*P<I*J5.KSl>  ♦03*P  ( I  » J5 ♦ *52 ) 

P  6  =  0 1  *  P  (  I  *  Jb  *  *S )  *D?*P(  I  t  Jb • K  S 1  )  ♦L>3*P(1*J6*KS2) 

IF  (  J.EU. JBl-i .OP.J.FU. OBJ >  GO  lo  181 

IF  ( J.EQ. JB2-1 .00. J.EO. JR?)  GO  TO  1182 

IF  ( J.Eu. JB3-1 .OP, J.EQ. J83)  60  10  182 

GO  TO  1185 

JK  =  J81 

PLC=PLCBl 

GO  To  183 

JK  =  JB  2 

PtC=PLCB2 

GO  TO  163 

JK=JB3 

PLC=PLC«  i 

COnT INUE 

PA=(?,*P( I »  JK ♦ 1  *  K  S  )  *  <  1. -PLC ) *(P<I*JK-1»KS) 

1  -2.*p < I «JK,K5) ) ) / n . ♦PIC)  ,  . 

PH= (?.*P < 1  * JK* 1 .KS1 )  ♦  ( 1  .-PLC) *  IP ( I  * JK-l .KS1 ) 
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c 


1  -2.*PU*JK,KSin  )/<l.*PLC> 

PC  = <  2»*P ( I ,  JK*1 »KS2) ♦ ( 1 .-PLC>«  <K< I , JK-1 »KS2) 

1  -2.«P<I,UK,KS?) > )/(l.»PLC> 

IF  (J.EQ, JK)  GO  TO  184 
P6=Dl*HA>02*P8^D3*PC 
GO  TO  1185 

184  P5=01*PA*U2*P8»D3*PC 

P4A=P  < I  * JK-1 »*S>  *3.* (PA-P ( I » JK»KS) ) 
P8B=Pa»JK-l»KSl)*3.*(P8-P(I*JK*KSi)  ) 

PCC=P( I , JK-1«KS2) ♦3.*(PC-P<I*JK*KS2) ) 
P6s01*PAA*D2*P8B*D3*PCC 
118*;  OPDS=Al*PZ*A2*P3*A3*P4 
0PD2sCl*PZ*C2*Pl ♦C3*P2 
U=CSA*CUl*0Pnb*CU2*HWri*0P0Z 

V=Cvl*OPO$*CV2*  <8l*P7*B2*P5*83*H6) ♦CAPGTB*HW0*DPDZ 
W=SNA*CW1*HWB«DP0Z 
IF  (J.EU.l)  V=0.0 
(IX*  1  ,-u*u-v*v-w*w 

CALCsA 1 2* ( ( 1 .♦0.20*QX/AI2) **3.5-1. ) /0.70 
IF  (MSURF.EU.2)  GO  TO  185 
ZXL ( I » J)  =  (W-V*DZ0YL ( I ♦ J) ) /U 
CPL(I»J)=CALC 
GO  TO  200 

185  ZXU<I»J)*<W-V»DZDYU(I*J> >/U 

CPU ( t » J ) =CALC 

200  CONTINUE 


1 

310 


JROOTX=JHOOT-1 
WRITE  (6,301)  CPSTAR 
00  310  J=1 , JROOTX 
YY=2,*Y ( J) /BSP AN 
WRITF  (b,302)  J , Y Y 
00  310  I=1N0S£,IMX 
TEMPsXCAPU)  >0.5 
UZDX1=UZ0XU(I,J) 

0Z0X2=0Z0XL ( I , J) 

WRITE  (6,303>  T£MP, CPU < I »U> »ZXU ( i *J) ,DZDX1 ,CPL ( I ,J) , 
ZXL(I.J) ,0ZDX2 

CONTINUE 

00  320  JsJROOT,JTIP 
YY=2.»Y ( J) /BSPAN 
WRITE  (6,302)  J»YY 
00  320  I  =  1LE » ITE 
TEMPsXCAP(I) ♦O.S 
r)ZDXl*DZDXU(I»J) 

0ZDX2=0Z0XL(I»U) 

WRITE  (6,303)  TEMP , CPU ( I , J) , ZXU ( I »U) , OZDX 1 »CPL ( I , J) » 


320  CONTINUE 


ZXL(I,J) ,DZDX2 


RETURN 

C 

301  FORMAT 
1 

30?  FORMAT 
1 
2 

303  FORMAT 


( IH-//////1X30HSURFACE  PRESSURE  DISTRIBUTIONS/// 
1X8HCPSTAR  =  E13.5///) 

( 1 M" ,17 HSPAN  STATION  J  *  I4,6X6H2r/8  *  F9,4// 
3X9H(X-XLE) /C, 12X3HCPU*8X7HSLOPE  U , 1 0X5HDZDXU , 1 2X3HCPL * 
8X7HSLOPE  L» 1 OX5HDZOXL/) 

(F12.4*2(F15.4,2F15.6)  ) 


END 
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